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I. INTRODUCTION 


A. BACKGROUND 


In 1906, ten years prior to the first appearance of the tank at the Battle of 
Cambrai, an Austrian munitions manufacturer created the first armor piercing round 
for use in small caliber rifles and machine guns [Ref. 1]. At the time, the only combat 
systems employing armor were warships, so most of George Roth's peers thought 
him crazy. But early in World War I, ingenious soldiers, seeking protection from 
bombardment and waves of assaults, began using steel plates for protection. In an 
attempt to penetrate these plates, the Germans used Roth's patented rounds. In 
1916, a new weapon, the tank, was employed on the battlefield and with its advent, 
munitions companies designed rounds that would defeat this new system. One of the 
first was created by Mauser. It was a 13mm rifle-fired projectile weighing 51.65g and 
had a muzzle velocity of 792 meters per second [Ref. 1]. 

With the passing of World War I, formal armor and mechanized units were 
created in the armies of the world's powers. The coming of age of the tank required an 
anti-tank weapon. In the 1920s and 1930s the general trend was to create weapons 
that fired larger projectiles or had greater muzzle velocity or both, exploiting the 
theory of kinetic energy. By the mid 1930s, the Germans and Poles each had weapons 
firing anti-tank projectiles in excess of 1,200 meters per second [Ref. 1]. 

Despite the heavy focus on high explosive (shape charge) munitions and the 
rocket launched projectile, the 1940s brought a significant breakthrough in armor 
piercing technology. The British formed a projectile using a dense tungsten core 
covered with a light alloy. The 57mm round was fired at 1,076 meters per second [Ref. 
1]. Additionally, the British Armaments Research and Development Establishment 
recalled research performed by the French using discarding sabots on armor piercing 
rounds. The British exploited this research creating a 20mm, 1.47kg round capable 


of attaining a 1,234 meter per second muzzle velocity [Ref. 1}. This projectile type 


became the primary armor piercing throughout the remainder of the war. 

Following World War II, anti-armor technology was primarily focused on anti- 
tank guided missiles. It was not until the 1960s that a new development in armor 
piercing technology occurred. The Soviets believed that the rifling required to main- 
tain aerodynamic stability usurped much of the energy of the propelling charge. Their 
solution was to use a smooth bore gun that propelled a fin stabilized dart. The dense 
sub-projectile was surrounded by sabots to fit the bore and capture the propellant 
gas behind it. This development reduced the bore friction allowing larger charges to 
be used and therefore achieve greater muzzle velocity. 

By the 1970s, every mechanized army used the armor-piercing fin stabilized 
sabot (APFSDS), even those with rifled bores. The tungsten core was replaced by 
depleted uranium. Depleted uranium offered a similar density to tungsten carbide 
and was a by-product of the nuclear industry. Since this period, small modifications 
have been made to the kinetic energy round, however it is, in essence, the same 
projectile. The majority of improvements have been made to the delivery system, 
particularly the fire control system. By having an in depth understanding of the 
ballistic characteristics of the projectile, the fire control system can compensate for 
varying parameters (atmospheric conditions, gun tube wear, cant, etc). 

Since the creation of the armor piercing round, many advancements in the 
shape, composition, and delivery have been made. However, one fact has remained 
unchanged; the kinetic energy round is an unguided, direct fire weapon. With “smart” 
munitions becoming increasingly common on the modern battlefield, the requirement 
for a new class of kinetic energy projectile has arisen. This new class of projectile 


may travel at very high speeds and require added maneuverability. 


B. STATEMENT OF THE PROBLEM 


The motivation behind this research is to attempt to increase the accuracy 


of armor piercing projectiles with minimal additional munitions manufacturing cost 


and with minimal modification to the launching platform. The main concept is to 
gain some of the accuracy associated with a “smart”, guided missile with a relatively 
inexpensive “competent” kinetic energy round. This concept may be built upon in 
the future to allow for “smart” kinetic energy weapons. 

Projectiles with unconventional shapes have been developed by British re- 
searchers. These projectiles are elliptic cones intended for travel at hypersonic speeds. 
In addition to the general elliptical shape. projectiles with added stabilizing surfaces 
have also been proposed, specifically, the addition of strakes and/or a tail flare. These 
shapes, designed for gun launched maneuver, are being analyzed by the United States. 
the United Kingdom, Canada, and Australia in accordance with a Key Technical Area 
(KTA) agreement. The United States Army Research Laboratory (ARL) is responsi- 
ble for the computational fluid dynamics (CFD) analysis of several of these projectiles 
in the preliminary design stage. 

In fulfillment of ARL’s responsibility, this study S analvze four proposed 
projectile configurations in steady flow of an ideal gas. Computations will be per- 
formed at Mach 5, 6, 7. and 8.2 at angles of attack of 0°, 5°, 10°, and 15°. The 
goal is to predict the aerodynamic characteristics of each projectile and to provide a 
comparative analvsis of the configurations. More specifically, the aerodynamic effects 
of the addition of stabilizing surfaces will be analyzed. The computational results 
will be compared to experimental results at Mach 8.2 and to computed data from 
previous research. Finally, this study will provide insight on the direction of future 


investigations into this class of weapon. 


С. PROJECTILE CONFIGURATIONS 


The four projectiles of interest in this study are based on a standard config- 
uration, a half power law elliptic cone measuring 540mm in length. This “baseline 
configuration is free of stabilizing surfaces. It is designated the H4 (H3 in previous 


studies [Ref. 10]) and its dimensions are depicted in Figure 1. The second design 
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igure 1. The H4 Projectile 
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adds four strakes to the H4 design. Each strake 1s tapered from the tail to zero at 
the nose following a constant circumferential angle. The strakes are symmetrically 
located on the body. This projectile is designated the H4A and is shown in Figure 2. 


The next configuration is identical to the H4 for the first 430mm in length. In the 
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Figure 2. The H4A Projectile 


110mm tail section, an elliptic flare is added. This flare is an elliptic conic section 
that increases the length of the semi-major axis at the base from 30mm as shown in 
Figure 1 to 35mm in Figure 3. The last configuration is a combination of the H4A 
and the H5, that is, a body with both strakes and a flare. The strake profile. identical 


to that of the H4A, blends into flare. Therefore. the strakes taper to zero (in height) 
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at both the nose and tail of this projectile. This configuration, the H5A, is shown in 


Figure 4. 





Figure 4. The H5A Projectile 


D. COMPUTATIONAL TOOLS 


The primary solver used in this research is the Zonal Euler Solver (ZEUS) 
created by the Naval Surface Warfare Center for tactical missile design. This solver 
uses the inviscid Euler equations and a spatial marching technique to numerically 
approximate the flow field about the projectile. From the flow field characteristics, 
the aerodynamic forces and moments can be determined. In steady, supersonic flow, 
disturbances are only propagated downstream. Therefore, given an initial flow field 
on a cross field plane near the nose, the remaining flow field can be predicted by 
advancing (or marching) to a cross field plane further down the projectile axis. This 
type of solver requires supersonic conditions everywhere, a condition easily satisfied 
under the desired parameters of this study. The theory of this solution scheme ıs 
discussed in detail in Chapter II. This solver has a reputation as both an efficient 
and robust program [Ref. 2. 3, 4]. Implementation of this code is facilitated by the 
creation of the ZEUS Graphical User Interface (GUI) by Edge and Clarke at ARL 
[Ref. 5]. 


To obtain the initial flow field on a cross flow plane mentioned above, a nose 
solver was used. The blunt noses of the projectiles precluded the use of the nose 
solver included in the ZEUS package. However, a ZEUS-coupled version of the Con- 
formal Mapping 3-Dimensional Transonic (CM3DT) code was used. This solver was 
originally designed for missile re-entry. It uses conformal mapping to produce a 
computational grid and solves the time-dependent, inviscid equations of motion in 
non-conservation form using a characteristic based finite difference scheme. By tak- 
ing the time asymptotic limit, a steady solution is obtained and used as the nose 
boundary condition in ZEUS. Both CM3DT and ZEUS rely heavily on slender body 
theory and the assumptions of inviscid flow. In particular, at hypersonic speeds, the 
entire area of interest lies between the body and the bow shock wave. Slender body 
theory allows the solvers to establish boundaries based on these conditions. Addi- 
tionally, the speeds at which the projectiles are being analyzed cause the boundary 
layer to become very thin, thereby, minimizing boundary layer effects. 

To facilitate analysis of the results, researchers at ARL have written scripts 
for use with NASA’s Flow Analysis Simulation Toolkit (FAST). These scripts enable 
the user to visualize the results obtained from both ZEUS and CM3DT. Additionally, 
prior to conducting computations, the scripts allow visualization of the computational 
grid, key in avoiding wasted computational time caused by poorly configured grids. 

All computations were performed on Silicon Graphics Octane workstations at 


the Naval Postgraduate School in Monterey, California. 


E. PREVIOUS RESEARCH 


The concept of non-circular munitions is far from novel. Research on bodies 
with elliptic cross-sections can be traced to the early 1950s. In Jack Nielsen's 1960 
landmark work on missile aerodynamics, he dedicates a section to the aerodynamic 
theory of these shapes and discusses their aerodynamic advantages [Ref. 6]. 


Specific work began on power-law bodies in the 1960s and the advantages 


of these shapes over circular cross-sections is well documented. Graves |Ref. 7| 
conducted a comparative study of the aerodynamic differences between a circular 
and elliptical (3:1 axes ratio) body. He concluded that the elliptical body generated 
approximately 25% more normal force than the circular body at supersonic speeds. 
This result was virtually independent of the angle of attack. Grave’s study also 
showed that the elliptic shape had greater lateral and directional stability which was 
accentuated at higher angles of incidence. 

Fournier et al. [Ref. 8] noted that increasing the ratio of the major axis to the 
minor axis increased lift. The research also showed that increasing Mach number (up 
to 4.63 in their research) had little effect on the Cy, — o curve. piss they stated 
that the Cy, — a curve was non-linear with the slope increasing with o. Kontis et al. 
[Ref. 9] maintained that this slope is further increased because of the separation of 
flow on the leeward side creating increased suction. 

As the concept of a new class of kinetic energy projectile came to fruition, 
research on elliptic projectiles was revisited. Recent published research has been 
conducted primarily by British academics and researchers from the United Kingdom’s 
Defense Research Agency (DRA). Of particular interest are the work by Edwards and 
Roper [Ref. 10] regarding the H4 and the works by Kontis et al. [Ref. 9, 11] comparing 
the H4 and H4A. 

The work by Edwards and Roper analyzes the H4 computationally from Mach 
9 to Mach 9 using a parabolized Navier-Stokes solver, CHASM, and the US Air Force’s 
Missile DATCOM (theoretical techniques and semi-empirical approximations). The 
computational results were validated using experimental data obtained at Mach 8.2. 
This data will be used in this work to validate the accuracy of the solutions obtained 
by ZEUS. 

The first work by Kontis et al. [Ref. 9] investigates the application of strakes 
to elliptic bodies. The H4 and the H4A are tested at Mach 8.2 at different angles 


of attack. In addition to the results of the gun tunnel experiment, schlieren pictures 


were used to determine the location of the flow separation at different angles of attack. 
The results obtained by Kontis et al. [Ref. 9] are used here for comparison. 

The other work by Kontis et al. [Ref. 11] analyzes a right elliptic cone with 
and without strakes in hypersonic flow. In this work, the authors conclude that in 
flow fields generated by these shapes, numerical studies can often produce valuable 
results despite being in poor agreement with experimental results. Their main con- 
tention being that although the numerical and experimental data did not agree, the 
computational data offers insights on the aerodynamic characteristics allowing better 
focus on further research topics. 

At the time of the publication of this research, no experimental or computed 
data had been published on shapes similar to the H5 or H5A. However, numer- 
ous works have been conducted on hypersonic cone-cylinder-flare configurations. Of 
particular interest is the work by Hughson and Chapman |Ref. 12] comparing the 
hypersonic aerodynamics of two flare tail bodies, one with a cone nose and the other 
with a blunt spherical nose. This research found that the difference in bow shock 
waves caused by the nose shape leads to very different flow conditions at the flare. 

Two additional studies pertaining to flare stabilization were applicable. The 
first by Guidos [Ref. 13] analyzes multiple cone-cylinder-flare configurations at Mach 
4.5 to Mach 9.0 and at small angles of attack. Comparisons of the effects of flare 
length and angles are made. Noteworthy is the drag-stability trade-off. He notes 
that flared projectiles incur a severe drag penalty to maintain the stability levels of 
a fin stabilized projectile. However, at the desired speeds, the fins of the current 
kinetic energy rounds are expected to fail because of excessive loads and heating. 
The other work by Danberg et al. [Ref. 14] also compares multiple cone-cylinder- 
flare configurations including a flare with strakes (no strakes on the forebody). In 
this work particular attention is given to the effects of flare design on drag, stability, 
and force and pressure distributions. 


The final work deserving note is that of Priolo and Wardlaw [Ref. 15], the co- 
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creators of ZEUS. In this study, three elliptic bodies were analyzed using experimental 
data and ZEUS. In this work, Priolo and Wardlaw describe techniques used in the 
application of ZEUS to bodies of elliptic cross-sections. The techniques used in Priolo 
and Wardlaw's work were used as a basis for the application of ZEUS in this study. 
Additionally, Priolo and Wardlaw document limitations of ZEUS when applied to 


non-circular cross-sections and describe the impact of these limitations on the results. 


F. SCOPE AND LIMITATIONS 
Consistent with the design methodology proposed by Moore [Ref. 2], this 


study is a preliminary analysis of potential projectile configurations. This research is 
meant to provide a basis for further investigation or elimination of untenable designs. 
Because of the preliminary nature of this research, an efficient, solver was utilized al- 
lowing a blend of robustness, ease of use, and accuracy. The use of this solver incurred 
certain limitations. The solutions achieved are for steady, uniform flow conditions. 
Performance analysis in complex flow fields requires in depth study beyond the scope 
of this work. Additionally, the Euler equations (used in the solvers) cause an inherent 
neglect of the effects of the boundary layer (skin, friction, aero-heating, etc.). Inves- 
tigation into these effects are left for future study. Finally, ideal gas conditions were 
used in this analysis, i.e. the effects of atmospheric conditions were neglected. 

Because of time limitations, only the effects of changing angle of attack were 
analyzed. 'This research may be used as a basis for further study into the effects of 
roll or yaw changes. 

Finally, because of the lack of experimental data for the H5 and H54A at the 
time of publication, the computed results are not validated. Based on the findings in 
this and other research reports, further experimental investigation may be conducted 


on these configurations. 
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THIS PAGE INTENTIONALLY LEFY BEANK 


П. MATHEMATICAL FOUNDATIONS OF 
THE SOLVERS 


A. CM3DT 
1. Background 


The CM3DT code was created in 1979 by Darryl H. Hall [Ref. 16]. In its 
original form, the CM3DT code was designed to predict the flow field over missile 
nose tips during atmospheric re-entry (altitudes less than 50, 000 feet). It was found 
that during re-entry, the nose tips ablated. This ablation, coupled with the supersonic 
and hypersonic speeds at which the object traveled, required a solver that could treat 
blunting, its associated shock waves, and any potential embedded shock waves. 

Under the conditions noted above (supersonic flow, low altitude, large Reynold’s 
number), the shock is almost completely inviscid except for a thin boundary layer 
adjacent to the body surface. By assuming inviscid flow, the aerodynamic character- 
istics can be predicted by the numerical integration of the time dependent inviscid 
equations of motion. Additionally, at low altitudes, the inviscid pressure distribution 
provides an adequate method of predicting the forces and moments (except in the 
axial direction). It was also determined that this method could be used to solve for 
the steady flow field. By integrating the time dependent equations over time, an 
approximation to the steady solution can be found. 

As inviscid Euler afterbody solvers were developed, the CM3DT code was 
implemented to solve the steady flow problem and provide steady initial conditions 
for the afterbody solvers. The code used in this research was modified specifically for 
use with the ZEUS code which does not have a built-in blunt nose solver. 

CM3DT takes several steps to obtain a solution to the time dependent prob- 
lem. First, CM3DT, through a series of conformal transformations in the projectiles’ 
meridional planes, creates a new coordinate system that is closely aligned with the 


body surface. The coordinate system facilitates the prediction of the bow shock lo- 
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cation and velocity and the formulation of the mesh required for the application of 
a numerical scheme. This transformation is particularly important for non-simple 
or slender nose geometries such as those in this study. When this step is accom- 
plished, a characteristic based scheme, the A—differencing method [Ref. 20], is used 
to integrate the transformed, time-dependent equations of motion. This numerical 
technique applies a modified version of the MacCormack finite difference scheme and 
achieves second order accuracy. 

When the time dependent equations are iteratively solved (in time), it is pos- 
sible to achieve a reasonable representation of the steady flow solution. This code 
finds the solution for steady flow by taking an asymptotic limit of the unsteady flow 
solution. 

This section describes the mathematical foundations used in CM3DT as ap- 
plicable to this research. This description will omit ideas not pertinent to this partic- 
ular problem. For a more detailed description, Hall [Ref. 16, 17] should be reviewed. 
The theory presented in this section is based on the research conducted by Moretti 


[Ref. 18, 20] and Hall [Ref. 17. 16, 19]. 


2. Governing Equations іп Cylindrical (r. y, 0) Space 
The equations used in this code are the 3- D, time-dependent, inviscid equations 
of motion and based on the coordinate system shown in Figure 5. The equations are 


applied in the non-conservation form 
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Figure 5. CM3DT Coordinate System 


where P is the natural logarithm of pressure, p, y is the isentropic exponent, s is 
entropy, and U, V, and W are velocity components in the z, y, and ó directions 
respectively. In these equations, the dependent variables are P, s, and the velocity 
components. The logarithm of pressure was chosen instead of pressure to minimize 
the large (orders of magnitude) changes in pressure. To close the system of equations. 


Hall uses the ideal form of the equation of state 


i = 


p= pie T” (11.2) 


where 
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and R is the gas constant. To achieve a solution to the steady flow problem, boundary 


бл 
8 


and initial conditions must be specified. Although any initial flow field will, in theory. 


converge to a unique steady state solution, to aid in convergence, Hall uses a technique 
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to provide a “reasonable” estimate. This technique is based on a user selected shock 
standoff distance which provides necessary shock information at the nose tip. In 
addition to the initial condition, two boundary conditions are required: the condition 
at the body surface; and the condition at the shock. The flow field on the body surface 
is determined numerically by applying the inviscid kinematic boundary condition (no 
velocity component normal to the body surface). The condition at the shock is 
found by numerically integrating the shock acceleration equation for the velocity and 
position of the wave. This information 1s unknown a prior? and 1s determined during 
the numerical integration of the problem. It should be noted that no downstream 
boundary condition is required based on the assumption that the flow across this 
boundary is supersonic. According to slender body theory, there will be no upstream 
influence and, therefore, no requirement for a downstream boundary condition. 

The equations and conditions described above are not applied in these forms. 
Hall uses a transformation technique developed by Moretti In 18] to create a more 
appropriate coordinate system. Consequently, these equations and conditions are 


mapped to a new system based on these conformal transformations. 


3. Transformation of the Coordinate System and Con- 
struction of the Computational Grid 


CM3DT applies a series of transformations such that the cylindrical coordi- 
nates (x,y, 0) are mapped to a new, near-cylindrical system (£,7,@) that is closely 
aligned with the body surface of the projectile. In this transformation, o 1s mapped 
directly to 0. Therefore. the true impact of this transformation occurs in the mapping 
of (z,y) to (£,9). In the transformed space, £ represents the direction aligned with 
the surface of the body surface while 7 is normal to body surface. This orthogonality 
is guaranteed by the conformality of the mapping of the orthogonal space (z,y). It 


should be noted that this transformation is time-independent. so t is mapped directly 
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to 7, a time coordinate. The transformations can be represented generally as 


€ — f(z.y. e) 

n = Q(z, y, 0) (11.4) 
0 — 

T —Í 


A representation of the untransformed and transformed coordinate systems (for any 
given $) is shown in Figure 6a. To accomplish the transformation in each plane, 
Hall uses Moretti's “hinge point” concept [Ref. 18] extended to three dimensions. 
In this method, a finite number (7) of hinge points are generated such that the 
points approximate the meridional cross-section of the object. These points simulate 
the body contour in the (z,y) plane, that is, in each @ plane. CM3DT generates 
equally spaced hinge points based on the wetted length of the projectile. The hinge 
points are a uniform distance (user input) from the body surface. These hinge points 
create corners (see Figure 6a). First, using a series of transformations, each corner 
is sequentially smoothed (starting from the nose and working aft) to create a nearly 
horizontal line in the (€,7) plane. Although each transformation is applied to all 
hinge points, its effect is to map the next off-centerline hinge point to the horizontal 
axis without affecting the upstream (previously mapped) points on the horizontal 
axis. This process is depicted graphically in Figure 6 [Ref. 16]. Before continuing, 


the following notation is defined. Each hinge point is described as 
hag (11.5) 


where the subscript denotes the i^ hinge point and the j'^ transformation. Addition- 
ally, each space is defined in a complex plane such that the original space is defined 
as 


21 = т + у (II.6) 


and 


zy 1 +Y (11.7) 
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Figure 6. CM3DT Hinge Transformations 


defines the space when the body surface is mapped to a horizontal line. J 1s the num- 
ber of intermediate transformations required to map all hinge points to the horizontal 
axis. The CM3DT code requires that the first hinge point, hı j, and the second hinge 
point, h2;, both lie on the projectile centerline. The first point must lie outside of the 
nose surface and the second lies inside the surface (see Figure 6a). Since these first 
two points are on the horizontal axis, they do not require transformation. However, 


the remaining points will require a mapping. Therefore, the number of intermediate 


mappings can be defined in terms of the number of hinge points 
J=1-2 (11.8) 


To accomplish the desired effect, the general form of each intermediate transformation 


is defined by 
== (2; = А15) —1 (II.9) 


where 
= 


ó; = 


4; — arctan 


mihj+25 — hj) m 
a о Ды 
Figure 6 depicts the complete transformation process. When all points are mapped 

to the horizontal axis, the space is zJ. 
Now that the body surface has been mapped to the horizontal axis, it is useful 


to have this mapping perpendicular to the image of the centerline (the straight line 


between hy, y and ha J). Hall accomplishes this with a simple square root transforma- 


Zj hs (11.11) 


One complication arising from the previous transformations must be resolved. Since 


tion 


the first two hinge points lie on the centerline, their location is independent of 0. 
The transformations thus far have been independent in each ¿ plane. Therefore, it is 
possible (and probable) that the images of the first two hinge points may be different 
for different ¢. To ensure that all mappings of these points coincide, the following 


stretching is performed 


Ç = G42J+1 a) 
where 
_ шм 11.13 
а = е (11.13) 
1J+1 


and hi ,,, represents the first hinge point in the J 4- 1 space and in the k^ ¢ plane. 


Having completed this series of conformal transformations, a new coordinate system 
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is defined where 
(= ш (11.14) 


and since 6 was directly transformed to б, 


EOE A (11.15) 


is complete. 

Knowing that the goal is to define a coordinate system that facilitates both the 
creation of a mesh and the application of a finite difference scheme, additional trans- 
formations are performed. Presently, the coordinate system (£,7, 8) is cylindrical-like. 
To create an equally spaced mesh, a rectangular system is required. Transformation 
to a Cartesian coordinate system also ensures easy numerical computation of deriva- 


tives. The transformation 
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where b(€, 0) defines the body surface points, c(£,0, 7) defines the bow shock position, 
and £r 1s the user defined downstream boundary. Body surface points are defined 
from a user input geometry file and the shock position is determined by numerical 
integration of the shock acceleration equation in time. This transformation has a 
normalizing effect on the coordinate system. That is, X becomes proportional to 0, 
Y varies between 0 (at the nose) and 1 (at the downstream boundary), and Z varies 
between 0 (on the body surface) and 1 (on the shock). 

50 far a series of intermediate conformal transformations have been performed. 


Since each of these mappings are conformal, the series of mappings results in a globally 
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conformal transformation. It should be noted that these transformations are time 
independent. Additionally, although the conformality of the global transformation 
guarantees that the orthogonal space (x,y) results in (€,7) being orthogonal, there 15 
no assurance that (€,7,0) is orthogonal. 

Having defined the computational coordinate system. the mesh can be con- 
structed. To build the grid, the user defines intervals in the X, Y, and Z directions. 
With these intervals established, the grid points are taken equally spaced in 7, €, and 


0. Figure 7 shows the grid generated for the H4/H5 nose at Mach 5 and a = 0°. 





Figure 7. Example of a CM3DT Generated Grid (H4/H5 at a = 0%, Mach 5) 


The change in the coordinate system affects the governing equations. In the 


next section, this effect is discussed. 





4. Transformation of the Governing Equations 
With a new coordinate system. it 1s necessary rewrite the governing equations. 
To convert from (7,y,0) {о (6, 7,0), the following chain rule relations are applied to 


(IL1) 


TEE 
22 DE a 
LAE 
д Sy Be ШЕ} 
(ШЕШЕ 
= am E 
до 7906! 90 
eO 
dt дт 
Now letting 
g= — = — = Gcos(—w) — iG sin(—w) (11.18) 


and using the fact that the transformations are conformal the Cauchy-Riemann con- 


ditions are applied 


fr = Ny = Re(g) (11.19) 
& = —n, — Im(g) 
Similarly, by defining 
ð lo | Ее 
= OE = ó, + ión (11.20) 
and transforming the velocity components 
и = U cos(w) — V sin(w) 
v = Usin(w) + V cos(w) (11.21) 
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the equations of motion and entropy can now be written as 
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To transform (11.22) to the computational space (.X,Y. Z), again the chain rule is 








applied 
9 — 9 +. 7 o 
ôr OT [OF 
д д д 
е == Z 5 
TI OY 07 
(II.24) 
д O 
B RI ae: 
On ОЙ 
O д 2, 
ex жЕ е: а ы н 
6х әү + к 
where 
Я 1 
Ag — D 
1 
у= 
` &(0) 
los —Y Yir, 
PEE (11.25) 
m 


Ze = —Z,|(1 — Z)b, + 2с] 
“ —Z,[(1 = 2)06 3r Z co| 
Z, = = ZZ 6; 


and band care the body and shock positions and their derivatives represent the body 


and shock slopes. These relations allow the governing equations to be written in the 


computational space 
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5. Applying the \—Differencing Scheme 

Having transformed the equations into a computational coordinate system, the 
A—differencing scheme can be applied. This method is a characteristic based scheme 
for hyperbolic problems in which finite difference approximations of the derivatives 
are taken in the directions determined by the sign of the characteristic slope. By 
using this type of scheme, the method ensures that the wave propagation phenom- 
ena is more accurately modeled. In mathematical terms, the use of A—differencing 
allows the domain of dependence information provided in the characteristic slopes to 
be incorporated into the finite difference scheme. To apply this method, the govern- 
ing equations are combined with characteristic compatibility conditions. ‘These new 
equations are integrated using a finite difference method and then the compatibility 


conditions are used to get the original form of the equations. Letting 
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the new form of (11.26) is 
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and the other coefficients are described by (11.28). X is the slope of the characteristic 
curves in the Z—T', Y —T', and X —T planes (subscripts denote the plane of reference). 
However, it is important to note that A has two values (denoted by subscript), both 
of which are positive for supersonic flow. If the flow is subsonic, Aı is negative and 
Aa is positive. The sign of A describes the direction of the domain of dependence at 
the point. For supersonic flow, only upwind information will be used in the finite 
difference derivative approximation of the spatial derivatives. However, for subsonic 
flow, both upwind and downwind information is used. A complete derivation of these 
equations and the compatibility conditions are described by Hall [Ref. 16] and Moretti 
[Ref. 20]. 

To integrate these equations, a modified version of the MacCormack method 
[Ref. 21] is used. This version of the MacCormack method is a one sided approach 
where in the predictor stage the backward and forward differences at nAz are respec- 


tively written 
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Similarly, the backward and forward corrected values are given by 
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where denotes the predicted value. Because of the corrector step, this technique 
achieves second order accuracy. For supersonic flow, only upwind information is 
passed, so backward differences are used. For subsonic flow, the terms related to 
the A, value use forward differences while the other terms use backward differences. 


The allowable time step size for this scheme is based on the Courant-Friedrichs-Lewy 


(CFL) criterion. For stability, the step size is selected such that the numerical domain 


of dependence is contained in the physical domain of dependence. 


6. The Steady Solution 

Using the methods from the previous section, the hyperbolic equations de- 
scribing the time dependent solution are solved. For steady flow, it is assumed that 
over some given time, the time dependent solution will approach the steady solution. 
Therefore, by applying a time-asymptotic relaxation, Hall uses the time-dependent 
equations to achieve a solution for the steady flow problem, a more efficient method of 
obtaining the steady solution. The time dependent equations are hyperbolic and can 
be solved by forward marching (in time) as a mixed initial-boundary value problem. 
Conversely, the time independent equations are elliptic in nature and must be solved 
as a boundary value problem. As the time dependent equations are integrated in 
time, an adequate approximation of the steady solution is achieved. To determine 
when this iterative process reaches a sufficient representation, Hall establishes four 
convergence criteria. When these four criteria are met, the CM3DT code ends con- 
sidering the final time iteration a close approximation to the steady solution. The 


criteria are: 
1. The stagnation pressure is essentially constant. It must fall within a predes- 
ignated tolerance of the theoretical value. 
2. The shock standoff distance must converge to a near constant value. 


3. The root mean square of the shock velocities must be decreasing and the 
magnitude must be within a predesignated tolerance of the magnitude of the 
free stream velocity vector. 


4. The total enthalpy at all points must be within à predesignated tolerance of 
the known steady state enthalpy (the steady state enthalpy is equal to the free 
stream total enthalpy). 


When each of these conditions are met, the steady flow field is solved and can be used 


in an afterbody solver. 


B. ZEUS 
1. | Background 


ZEUS is an Euler solver created by the Naval Surface Warfare Center to 
perform aerodynamic predictions on tactical missiles with body-wing configurations. 
oince its implementation in 1986, ZEUS has been used extensively on both axisym- 
metric and non-axisymmetric configurations. The efficiency, robustness, and relative 
ease of use compared to Navier-Stokes solvers has made it popular with project man- 
agers particularly in preliminary design stages [Ref. 2]. For these reasons and its 
ability to handle thin, low aspect ratio lifting surfaces, the Army Research Labora- 
tory has used it for aerodynamic predictions of projectiles. 

ZEUS uses a second order Godunov finite volume scheme to explicitly integrate 
the inviscid Euler equations. Solutions are obtained by marching spatially from nose 
to tail. Like CM3DT, ZEUS generates a grid using the body surface and bow shock 
as boundaries. However, the grid accuracy is enhanced by the ability to divide the 
computational domain into one or more zones. Each zone can be mapped separately 
based on the cross plane geometry. 

ZEUS requires several assumptions. First, ZEUS is a supersonic flow solver. 
The flow must be at least supersonic everywhere on the projectile such that distur- 
bances do not affect upstream flow. Additionally, ZEUS assumes steady flow and 
marches spatially solving for a steady solution. Finally, since the inviscid Euler equa- 
tions are used, the boundary layer is assumed to be very thin and its effects are 
neglected. 

In this section, an overview of the methodology used by ZEUS will be pre- 
sented. For an in-depth review of the theory, the reader should consult the work by 


Wardlaw et al. [Ref. 22, 23, 24]. 
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2. Control Volume Notation and the ZEUS Coordinate 
System 


Although mesh generation will be described later, it 1s important to establish 
the notation defining each control volume. In finite volume schemes, properties of 
cell edges (versus points) are calculated. Each cell is described by six edges. The size 
and location of each edge is defined by the location of a grid point which becomes 
the center of each cell face. The notation used in ZEUS for each point is (n.m, k) 
which is analogous to the ZEUS coordinate system (z,y,z). In ZEUS, z is up, y is 
left. and z points aft. Figure 8 [Ref. 24] shows the ZEUS coordinate system and the 


control volume notation. Having defined the cell notation, the Euler equations can 


Форте жы) ж 





Figure 8. Control Volume Notation and ZEUS Coordinate System 
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be defined and converted to a usable form. 


3. | Governing Equations 


The inviscid steady flow Euler equations for perfect gases are defined as 
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These equations, based on the conservation of mass, momentum, and energy, can be 


placed in integral form 
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In these equations, AŽ „ is the area of the cell edge containing (n, m, k), n = (nz, ny, n=) 
is a vector normal to the cell edge (see Figure 8) and has a magnitude equal to the edge 
area, and V = n,u+n,v+n,w. (u,v, w) are velocity vectors in the (x,y,z) directions 
respectively. U represents the flux passing through the shaded cell sides in Figure 8. 
That is, U is the flux passing through each z —constant cell face. Conversely, the F 
values represent the flux through the non-shaded sides. The systems found in (11.34) 
and (11.36) are in conservation form and are closed under the constraint based on the 
perfect gas equation of state and constant total enthalpy 
Pri 


1 
Ho = 一 +u tv tw’) (11.38) 
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Both equations (11.34) and (11.36) are valid for flows producing strong shocks. 





4. Grid Generation 
Having established the governing equations, it is now logical to define the 
computational grid. The grid generation in ZEUS differs from that of CM3DT 1n two 
Ways: 
e Onesimple transformation is used versus a series of conformal transformations 


e Grid points serve as defining points for control volumes 


ZEUS avoids the need for conformal mappings by allowing the user to divide each 


cross flow plane into quadrilateral zones thereby providing a method for dealing with 
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complex geometries. Typically, this zoning 1s used for wing-body configurations. In 
this work. all grids consist of a single zone. The zone boundaries in these cases are 
the body surfaces, the bow shocks, and the leeside and windside meridians. Figure 9 


[Ref. 24] defines the ZEUS notation for each zone. 


LEEWARO 
LEEWARD 
PLANE 
TO? FIN 
SURFACE 





SHOCK 





WINDWARD 
WINDWAHD 
PLANE 





BOTTOM 
FIN SURFACE 





Xa, Ya 


Figure 9. ZEUS Zone and Mesh Notation 


34 


ZEUS creates the mesh using the mapping T; - T» where 


BE (AO) = (Gre V) 


(11.39) 
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In Cartesian coordinates, this transformation becomes 
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where the points ø, b, c, and v are shown in Figure 9 and are described by 
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апа (s;, 7;) are coordinates of the cell corners. The functions f(£) and g(r) determine 


the clustering of the mesh. The user defines these functions by designating arc length 


percentages at each zone edge. For a uniform mesh, 


HE) =€ E 
T) 


In this study, the grid was clustered only in the circumferential direction. 


5. Transforming the Governing Equations 
Having established a computational grid based on the (£,7,¢) coordinate sys- 


tem, a final transformation of the governing equations can be performed. 


In order to evaluate the terms in (11.36), this system must be put in non- 


conservation form 
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It should be noted that Z is mapped directly to C. Derivation of this form of the 
Euler equations is performed by using a non-linear transformation of variables. The 


variables in this system can be evaluated using the numerical scheme described in the 


following section. 


6. The Numerical Scheme: An Application of a Second 
Order Godunov Method 


To numerically solve the non-linear systems,(II.36) and (II.44), a second order 
Godunov predictor-corrector scheme is used. This method determines the fluxes, F’, 
given in (11.36). This is accomplished by evaluating (11.44) using a twelve step method 
described by Wardlaw and Davis [Ref. 22]. For a more detailed (and theoretical) 


description of the method, Davis [Ref. 28] should be reviewed. 
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The first step in this process computes the partial derivatives of the primitive 
variables, (p, p, u, v, and w), with respect to the computational coordinates, (£ and 
1). These values are obtained using a limiter and by computing the central differences. 
The role of the limiter will be discussed in Chapter V. However, it should be noted 
that the value of the limiter is usually one for all interior points. Letting f represent 
any of the primitive variables and o represent the computational coordinates. an 


approximation of the partial derivatives is obtained using the equation 
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and n and m represent the n^ and m^ points in the y and € directions respectively. 


This computation is performed for each point in the z =constant plane. 
Similarly, the metrics, z, Ey, E2, Тг Ny, Nz, Gz, Gy, and ¢, are evaluated at each 
point. ‘These values are used in the predictor step and in determining the step size, 


Az. The values of the metrics are obtained analytically using the equations 
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and 


Details of this derivation are in Wardlaw and Davis [Ref. 
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22]. Having computed the 


partials of the primitive variables and the metrics, the step size can be calculated. 


where 


Ihe step size, Az, ıs determined by the CFL condition. 
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and a 1s the speed of sound. This condition ensures that the domain of dependence 


of the numerical scheme contains the domain of dependence of the Euler equations. 


ZEUS allows for user modification of the step size using a CFL safety factor. Appli- 


cation of this factor will be discussed later. 


At this point, the step size, the metrics, and the primitive variable partial 


derivatives have been calculated. With this information, the solution can be advanced 


and the predictor values of the primitive variables can be calculated. The calculation 


of the predictor values is taken at z + ==. The predictor values of the primitive 


variables, p, p. u. v. and tv. are calculated by substituting the metrics and derivatives 


into (11.44) and solving the system. 
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Before solving for the corrector values, information describing the control vol- 
umes and cells must be computed. The goal of this step is to calculate the coordinates 
of the control volume corners and to determine the vectors normal to the cell edges 
(see Figure 8). To determine the coordinates for the corners of the control volume, 
the grid generation transformation is used (11.41). The normal vectors for each cell 
edge are calculated using 
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where V, 4 is the vector from (2p, Yn, Zn) tO (Lm, Ym; Zm) and the subscripts correlate 
to the corners of each cell edge in Figure 8. 

Two critical pieces of the problem are known: the predictor values of the prim- 
itive variables and the location of the cell corners. From these pieces of information. 
the predictor values at the cell edges can be interpolated. However, each interior 
edge is influenced by two points. For example, in Figure 8, the lower edge of the 
cell defined by (n,m, k) is influence by both (n,m,k) and (n — 1,m,k). This will 
be resolved in the next step. To interpolate the values at the cell edges (using the 


notation in II.46), the following equations are used 
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To determine, the actual value of each property at the cell edges, the Riemann problem 
is constructed. For supersonic flow, the Riemann problem represents the confluence 
of two supersonic streams. At the intersection of these streams, both are turned 
to a common direction and have identical pressures. Figure 10 depicts pressure 
direction curves for the confluence of two streams. Solving the Riemann problem 
results in constant properties for both streams. By using this technique and solving 
the Riemann problem, the properties at the cell edges can be found. ZEUS has the 


ability to solve the complete Riemann Problem or to obtain an approximate solution. 
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Figure 10. Pressure-Direction Curves Depicting the Riemann Problem 


In this study, the complete solution was used. The setup of the approximate Riemann 
problem is thoroughly explained in [Ref. 22, 28]. 

The Riemann problem is solved by computing the final flow direction i.e. the 
direction that generates identical final pressures in both streams. The initial flow 
directions, ó- and ó+}, and the initial pressures, p- and p+, are related to the final 
direction and pressure, ó; and p;, by the shock and expansion equations given below. 
The first equation, the shock equation, is valid for py > pz, and the second equation, 


the expansion equation, holds for pz > ру. 
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where 





(11.56) 
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and ó+ and p+ represent one of the flow directions and pressures. To solve the shock 


reram MA 1 





and expansion equations for both flow directions, an iterative process is used. For 
the n+ 1 iteration, two points are known on each of the curves shown in Figure 
10. For the - stream curve, the points (67, p*) and (6?-!, p*-!) are known. A line 
is drawn through these points on each curve. The intersection of this line serves as 
the estimated pressure and is substituted into (IL.55) as pr. Similarly, p^ and ё" аге 
substituted for pz and ó+. Making these substitutions allows for the evaluation of 
0;.. Substituting the values from the + stream curve results in ôf}. The iteration is 
terminated if |óg, — 6;5..| « 1078. 
Before using the solution from the Riemann problem, the angular orientation 
of the cell edge must be computed. This direction is given by 
n 


®n+4,m 


mx 1 Tn, 


n+ 5 


0 = arctan — (TES) 


n+3,m 
where the subscripts on n describe the direction of the normal vector from the cell 
face. 

With the cell edge orientation known, the properties on the cell edge can be 
calculated. The pressure and density are found by using the Riemann solution in the 
@ direction. The magnitude of the velocity is computed using the results from the 


Riemann solution 
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Converting this value to a dimensional value, gg. the velocity components can be 


determined f 
_ qasin(ó)nx — u,ny 
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Ihe cell center from which v, is determined is based on the relative sizes of 0 and бу. 
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Having the properties at the cell edges. the corrector values are determined. The cell 
edge properties are substituted into the Euler equation in conservation form (11.36) 


and the fluxes are evaluated. The flux evaluations are decoded giving the final flow 


properties at the cell side at =z + Az. 
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where À, ,, is the cell face area defined in Section 3 and 
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= 
and u; corresponds to the i'^ row of U} „, (see (11.37)). 


7. | Boundary Conditions 


(11.63) 


Discussion of boundary conditions for this scheme has thus far been conspicu- 


ously absent. The boundaries in ZEUS are identical to those in CM3DT, specifically, 


the body surface and the bow shock wave. However, since ZEUS uses a finite volume 


scheme, no grid points lie on the boundary; cell edges lie along the boundary. The 


properties along these edges are based on extrapolation from the center point. It 


should be noted that a modified version of Godunov’s technique is applied to edges 


adjacent to the body surface. This scheme uses one sided differences to approximate 


derivatives normal to the surface. Additionally the flux passing through cells adjacent 


to the body is constrained by the tangent flow boundary condition. A full description 


of this modified version is in [Ref. 22]. 
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ПІ. MODELING THE PROJECTILES 


A. BACKGROUND 


As with any simulation, the accuracy of the results is based on the accuracy 
of the input. With this in mind, it was imperative to accurately portray the shapes 
of the four projectiles. Both CM3DT and ZEUS accept user generated geometry 
files in Plot3D (P3D) format. In this format, grid points aligned on axial (cross- 
sectional) and circumferential planes describe the surface and each point is stored as 
an (x,y,z) coordinate. MATLAB computer codes, h4generator.m, h4agenerator.m, 
hógenerator.m, hàagenerator.m, were designed to construct each projectile in P3D 
format. In each code, the points were output to a ASCII file in Plot3D readable 
form. The general format for loading points into this form is 

WRITE((X (i,j) ,i=1..imax),j=1..jmax) 
WRITE((Y(i,j),i=1..1max),j=1..jmax) 
WRITE((Z(i,j),i=1..1max),j=1..jmax) 
where the z represents the axial direction and 7 represents the circumferential direction 
as shown in Figure 11. These generation codes are located in Appendix A. 

The H4 was easily defined mathematically and so constructing an accurate 
model was simple. Similarly, the H5, being a smooth surface, was also relatively 
easy to construct. The projectiles with strakes offered the greatest challenge. First. 
modifications to the shape had to be made to make it usable in ZEUS. Then, a 
mathematical method for describing these shapes had to be developed. This chapter 
will discuss the construction of the geometry files defining the four projectiles. It 
should be noted that because of pitch plane symmetry associated with pitch only 


flight, the models used in these tests consisted of only half of the projectile (positive 
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Figure 11. Notation of Indices for Creation of Plot3D Files 


В. НА 


The H4 is easily defined in ZEUS Cartesian coordinates (see Figure 11) as 


232° ү Оу =1l5z 0<2 < 540 


(III. 1) 


where z = 0 is at the nose and z = 540 is the value at the tail. By applying the 


standard definition of an ellipse in parametric form 
у = асоѕ ф 
т = bsinó (op s 2m 
then (III.1) becomes 
250? sin $* + да? соѕ ф° = 152 
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To determine the magnitudes of the semi-major and semi-minor axes, b and a, re- 


spectively, the conditions ó — 0 and ф = > are applied. At ф = 0 
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(111-4) 


(111.5) 





Using a direct mapping of z to (, the H4 is parametrically defined as 


7 = a cos ó 
TS O (11.6) 
z = 


where a and b are defined above and 
(HI.7) 


From these equations, the model was generated using a MATLAB code. This algo- 
rithm modeled the H4 (half shell) with 60 circumferential points (z and y) coordinates 
and 109 axial (z) coordinates for a total of 6,540 points. Since the H4 serves as the 
basic shape for each of the other projectiles, the derivations noted above will prove 


crucial in modeling the remaining objects. 


C Но 

Like the H4. the H5 was smooth and relatively easy to define mathematically. 
For 0 € ¢ € 430, the H4 and H5 are identical. Therefore, this section will discuss the 
modeling of the flared tail that distinguishes the H5. 

The flare of the H5 was designed as an elliptical conic section. The ellipse 
that describes the tail of the flare has a different eccentricitv than that of the H4. 
Therefore. the flare was not easily defined mathematically. But because the base of 
the flare was described by a known equation, then for every о, the (x,y) coordinates 
could be determined for both the z = 430 and z = 540 cross flow planes. At z = 430, 
the beginning of the flare, these points were computed using (III.6). The tail of the 


flare was an ellipse defined by 


22.85?y* -- 35? x7 -— 799.757 (111.8) 
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Using the same method used in the previous section, the semi-minor and semi-major 


axes were found 


199.75 
а= 2978 = 39 
(111.9) 
799.75 
b = —— = 22.85 
35 


Having determined these constants, the parametric equation for the tail plane became 


у = 35cosó 
z=18csd O0< ¢6< 27 


(111.10) 


For every specified ¢, points (x, y,z)z=430 and (z,y,z)z=549 were determined using 
(III.G) for z = 430mm and (111.10) for z = 540mm. From the pair of coordinates 
(х,у, 2),-430 апа (=, у, 2), 540 іп а given ¢ plane, the equation of a line was derived 
and equally spaced grid points along this line were computed. This process was 
performed for each circumferential plane. Like the H4, the modeled H5 had 109 axial 
grid points (of which 21 defined the flare) and 60 circumferential points. 


D. H4A 


The projectiles with strakes offered two problems not found in the projectiles 
without strakes. First, the bodies were no longer smooth and could not be simply 
defined by mathematical equations. Even using a grid generator, the design of the 
H4A and H5A was cumbersome. Second, the angles between strake sides and the 
body were potential problems for the Euler solver, particularly on the shoulder side 
of each strake. ZEUS, like other CFD codes, relies on the orthogonality of the grid to 
obtain an accurate solution. In fact, ZEUS was unable to properly “march” axially; 
negative pressure values occurred in these areas. To overcome this problem, the shape 
of the strakes had to be modified while maintaining an accurate depiction of the true 
shape of the projectiles. A "smoothing" of each strake's shoulder side was performed. 


This modification was minor and resulted in only a 2% increase in the cross-sectional 
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area at the tail. Figure 12 shows the difference between the actual and modeled strake 


in a cross flow plane. 
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Figure 12. Actual and Modeled Strake 


In keeping with the methods used to model the H4 and H5, the H4A was also 
modeled circumferentially. Specifically, the projectile was modeled in five pieces: the 
body from ¢ = 0° to the upper strake (¢ = 0° in ZEUS coordinates is the z-axis); the 
strake; the body between the upper and lower strakes (the shoulder of the projectile); 
the lower strake; and the body from the lower strake to o — —180?. Because the 
strakes followed the ¢ planes, this was the easiest method of modeling the projectile. 
The body points were still easily defined by (111.6). The strakes were constructed by 
defining each side: the leeward or windward side; the top; and the shoulder side. 

The leeward and windward sides were modeled exactly, that 15, no smoothing 
was used. The angle between these strake sides and the body allowed a grid to be 


generated and ZEUS was able to march across these areas. 
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The top of the strakes were slightly curved instead of being squared off as 
designed. The equation for the flare tail, (III.10), was used. The reason for this will 
be discussed in the next section. 

Figure 12 shows the sharp edges caused by the strakes. The strake sides closest 
to the shoulder of the projectile have significantly steeper slopes than the leeward and 
windward sides of the strakes. These steep slopes result in smaller angles between 
the strake and the body. In an attempt to alleviate this problem, the strake sides 
were smoothed. A comparison between the actual and modeled strake side is depicted 
in Figure 12. The coordinates show the difference between the actual and modeled 
strake-body intersection points. 

In two dimensions (the cross flow plane), the strake sides are straight lines 
as shown in Figure 12. Because the strakes taper, the length and slope of these 
lines change in each cross flow plane. That is, moving toward the nose, these lines 
become shorter and have decreasing (in absolute magnitude) slopes. To smooth the 
entire side, the lines in each cross flow plane were represented as circular arcs. The arc 
sections in each cross flow plane were defined by a circle whose center and radius were 
dependent upon the axial station. To implement this method, three points defining 
each circle were required. The endpoints of each arc, A and B in Figure 13, were 
used. It is important to note that B, the point where the strake and body intersect, 
is not the exact location of intersection; point C describes where the actual strake 
intersects the body. Point B was chosen by selecting an arbitrary distance from point 
C, approximately 4mm. To compute a third point, thereby uniquely defining the 
parent circle of the arc, an approximation based on midpoints was used. 

The midpoints of segments AB and AC, denoted E and D respectively, were 
found. The midpoint of the segment connecting these points, F, was the third point. 


Now an equation for the circle defined by the points, A,B, and F, was found by 
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solving 
r= {т lua — b) 
r = /(rg —h)? + (yp — А)? (111.11) 
r= yer =h} + (ur — 0) 
where r is the circle radius and (h,k) is the center. Once the equation of the circle 
was derived, an arbitrary number of grid points were computed to represent the arc. 


This process was repeated for each cross flow plane. 





Figure 13. Construction of the Modified Strake Side 


As with the previously modeled projectiles, 109 axial coordinates were used. 
However, circumferentially, 61 coordinates were used. Additionally, the circumferen- 
tial planes were not uniformly distributed. To capture the smoothing on the shoulder 


sides and tops of the strakes, points were clustered in these areas. 


E. H5A 


The construction of the H5A model was a combination of the processes of 


modeling the H5 and the H4A. However, when the flare begins at z = 430, the 


ol 





| 
i 
| 


strakes begin to blend into the body. But as previously mentioned, the strakes follow 
a constant $. Therefore, the methods used to compute the sides of the strakes were 
identical to those used in the development of the H4A with different strake-body 
intersection points. These points were determined using the flare building algorithm 
of the H5. 

It was mentioned previously that the strake top was slightly curved (virtually 
impossible to detect visually) instead of the squared top in the design. Because the 
strake is supposed to blend completely into the flare at z = 540, the curved top made 
this transition easy to design with negligible aerodynamic effects. 


Figure 14 depicts the modeled projectiles. 





Figure 14. Models of the H4 (top left), H5 (top right), H4A (bottom left), H5A 
(bottom right) 





F. PROJECTILE NOMENCLATURE 

Having defined and named the projectiles, it should be noted that previous 
works have used other nomenclature to define these projectiles. For clarity, in this 
work, projectiles with strakes are indicated by the “A” identifier and application of 
the flare is denoted by the numerical identifier. ‘Table I shows the nomenclature used 


In previous works. 


Present Nomenclature | Prior Nomenclature 


H4 
H4A 
H5 not named 


H5A H5 


Elliptic Projectile 
Projectile with strakes only 
Projectile with flare only 
Projectile with strakes and flare 























Table I. Present and Prior Projectile Nomenclatures 
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IV. APPLYING THE CM3DT CODE 


A. BACKGROUND 

Since the H5/H5A differ from the H4/H4A only in the tail flare, nose solutions 
were required for only two geometries: the projectile with strakes and the projectile 
without strakes. That is, the nose solutions for the H4 and H4A were used for 
the H5 and H5A respectively. Since the method of obtaining the nose solutions was 
nearly identical for both geometries, the discussion applies to both nose configurations 
except where noted. To avoid confusion, the ZEUS coordinate system will be used 
when describing objects in this chapter. 

The application of the CM3DT code is actually a three step process involving 
manipulation of the model geometry files, the use of a geometry pre-processor, and 
then the running of the solver. 

In this chapter, the methodology, user inputs, and code output pertaining 
to the nose tip solution will be discussed. As noted previously, only information 
pertinent to this specific problem will be addressed. Additional information on the 
application of the CM3DT code to other problems can be found in the CM3DT user 
guides [Ref. 25, 26]. 


B. PREPARATION OF THE GEOMETRY INPUT FILES 

To prepare the geometry file to be input into CM3DT code, four operations had 
to be performed. First, since the CM3DT is concerned only with nose solutions, the 
nose of each projectile had to be parsed from the afterbody. Unlike many conventional 
(cone-cylinder) projectile shapes, the noses of the projectiles in this study were not 
clearly defined. In this research, the length of the nose was 5mm. To ensure all 
required geometric data was obtained, the first 10mm of the projectiles was used. 
The size was based on the desired location (output) of the solution. This will be 


discussed in further detail later in this chapter. 
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Next, the number of points defining the nose was rescaled. CM3DT limits the 
number of point to 162. Specifically, only nine circumferential and 18 cross-sectional 
planes can be used. The models described in the previous chapter inadequately define 
the nose cross-sectionally (only three planes instead of the 18 in the first 10mm) and 
over-define in the circumferential direction (over 50 planes). Additional problems 
existed. 

The version of CM3DT used in this study was modified specifically for use 
with ZEUS. However, one major idiosyncrasy remains; the ZEUS and CM3DT coor- 
dinate systems are different. The orientation of the projectile had to be modified to 
align with the CM3DT coordinate system. Both computer programs read standard 
P3D geometry files that provide grid points in Cartesian coordinates, however the 
orientation of the axes is different. The ZEUS coordinate system defines x as up, y to 
the right, and z running axially aft. Conversely, CM3DT defines y as up, z to the left, 
and z running axially aft. Figure 15 compares the two different coordinate systems. 
Since the majority of the calculations on these projectiles occurred in ZEUS, they 
were modeled in ZEUS coordinates as discussed in Chapter III. 

Finally, CM3DT failed (specifically the pre-processor) when the nose tip was 
the origin. The origin caused a singularity and the nose tip had to be modified. The 
singularity occurs when transforming from Cartesian coordinates to the cylindrical- 
like system described in Chapter II; this transformation requires the evaluation of =. 
A simplistic modification was made i.e. eliminate the origin. The nose tip of each 
projectile was “shaved” resulting in the nose tip being described by an extremely 
small plane instead of a point, the origin. In ZEUS coordinates, the nose tip is a 
plane at z — 10^? instead of the origin. The method for overcoming this problem 
was provided by Edge [Ref. 30]. 

Although any of these problems left unmodified would cause catastrophic re- 
sults, the required changes were easily implemented. The h4generator.m and ћ4 адепет- 


ator.m codes were modified to create two nose generating codes h4nosegenerator.m 
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Figure 15. ZEUS (a) and CM3DT (b) Coordinate Systems 
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and h4anosegenerator.m. 'The modifications were simple and accomplished four goals: 
to parse the nose from the projectile; to redefine the number of planes (grid points) 
describing each nose; to map the objects into the CM3DT coordinate system; and to 
remove the singularity caused by the origin. 

To parse the nose and redefine the grid points describing the nose, each z 
vector was changed from a vector of length 109 going from 0 to 540mm to a vector of 
length 18 going from 1071? to 10mm. Similarly, the number of circumferential planes 
was reduced to nine. To transform each point into the CM3DT coordinate system, 
the output of the generators was modified. In their original (ZEUS) form, the code 
outputs grid points as shown in Chapter III. 'To transform into CM3DT coordinates, 
the data was written to the geometry file in this format 

WRITE((Y (i,j) ,i=1..imax),j=1..jmax) 

WRITE((-Z(1,j),1—1..imax),j—1..jmax) 

WRITE((X(i,j),i=1..imax),j=1..jmax) 
where the indices are shown in Figure 11. 

Special treatment was used on the nose with strakes. Because of the limited 
number of circumferential planes (or points in each cross-sectional plane), the nose 
with strakes could not be accurately modeled; at least eleven points are needed. 
Using the nine points available, an approximation was constructed. In this model, 


the strakes are pointed instead of squared (see Figure 16). 


C. GEOMETRY PRE-PROCESSOR 

Once the modifications to the geometry were completed, a pre-processing code, 
preproc.z was run. This code, part of the CM3DT software package, simply converts 
the P3D file to a readable form for CM3DT. The Cartesian coordinates are converted 
to the cylindrical-like system (zx, y, $) described in Chapter II. When this code was 
run, the inputs in Table II were used in the required input file preproc.inp. It should 


be noted that the “1” usually found on the first line of a “multi-grid” formatted P3D 
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Figure 16. Actual (Solid) and Approximated (Dotted) Straked Nose Cross-Section 


file (the format accepted by ZEUS) must be removed prior to execution of the code. 
CM3DT uses the “single-grid” P3D format which omits thıs line. In the codes shown 


in Appendix A, the “1” ıs written to the output files. Visualization programs require 


this number in P3D files. 


The output from the preprocessor is cm3dt.geo, a file that 1s read during the 


execution of the next step of the process. 


isa Pitch plane symmetry (1=yes) ] 
Imax Number of circumferential planes | 9 
fname | Name of the file to be processed | h/nose.grd 


















Table II. Geometry Pre-Processor Input Variables 
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D. THE FLOW FIELD SOLVER 
1. Required Inputs 


The solver, Cm3dt2.x, requires two input files: cm3dt.geo (pre-processor out- 
put) and CM3DT.INP. The CM3DT.INP file consists of the flow conditions and other 
information described in Chapter IJ. The variables and values used in this research 
are shown in Table III. The variables below the solid line are variables unique to the 
ZEUS-coupled version of this solver. An example of CM8DT.INP and preproc.inp are 


in Appendix B. At this point, the choice of some of these values deserves explanation. 


Mach number 
Angle of attack 

Angle of sideslip 

Use of A-differencing scheme (1—yes) 

Number of time steps 

Number of grid points between surface and shock 
Number of grid points along each circumferential plane 
Number of circumferential planes 

Number of time steps between complete output 
Effective nose radius 

Location of nose tip in z direction (ZEUS coordinates) 
Points in z direction (ZEUS coordinates) describing 
downstream boundary 
Number of cells between edges 1 and 3 in ZEUS 
Number of cells between edges 2 and 4 in ZEUS 
Axial location where starting data is required 
Location in ZEUS coordinates of X9 in CM3DT coordinates 
Location in CM3DT coordinates of Z in ZEUS coordinates 



























Table III. CM3DT Input Variables 


The number of time steps, KMAX, was chosen based on a recommendation 
from Edge [Ref. 29]. This value had to be sufficiently large so that the convergence 
criteria described in Chapter II could be met. 

The density of the solution grid, described by NMAX, MMAX, and LMAX, 


was maximized to obtain the highest resolution possible. 
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The effective nose radius, Rw, is used to estimate the shock stand-off distance 
and the location of the hinge points. For blunt objects, the bow shock is separated 
from the object. This shock stand-off distance is based on the bluntness of the object. 
The default value of this variable was used. The value selected represents a reasonable 
stand-off distance and allows the hinge points to accurately simulate the contour of 
the body surface. The selection of the value of Ry is arbitrary, because of the iterative 
nature of this solver. In theory, convergence will occur for any value of Ry. However, 
a reasonable estimate speeds the solution process. 

XBB describes the downstream boundary i.e. the last axial plane for which 
CM3DT will provide a solution. Two requirements are imposed on this value. It 
must be less than the length of the nose and it must be larger than the location of 
the desired starting data. Because the solution is initially aligned with the body and 
therefore not a z =constant plane, the solution at this boundary cannot be used for 
a ZEUS initial flow condition. A built-in interpolation algorithm is used to obtain 
a solution field for a z =constant plane. A tenable value of XBB is determined by 
the shape of the projectile and the desired location for the ZEUS starting data. It 
is also important to note that CM3DT is accurate for nose solutions and not the 
afterbody; its utility in this study is solely to provide initial data for the afterbody 
solver. Therefore, the value of XBB was selected to accurately represent the shape 
of the nose yet allow ZEUS to provide the flow field solution for the majority of 
the projectile. Additionally, the value of XBB had to be sufficiently large to allow 
interpolation of flow conditions on the z =ZSTART plane. 

Of the inputs required for the ZEUS coupled version of CM3DT, NA, MA, 
and ZSTART are of interest. NA and MA describe the density of the output grid. 
Because of a conversion process used in ZEUS that will be discussed in the following 
chapter, any reasonable choice for these values can be used. As mentioned above, the 
value of ZSTART is chosen to obtain a sufficient representation of the projectile’s 


nose while allowing ZEUS to solve for the afterbody flow field. Since these projectiles 
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do not have a clearly delineated nose, the choice of this value was arbitrary. A decision 
was made to keep the length of the nose below 196 of the projectile length. Since the 
output at this value is only serving as initial solution for the marching solver, the 


location of the solution was approximated at one percent of the projectile length. 


2. Output 

The execution of Cm3dt2.x results in a text output file, (CM3DT.OUT), a 
ZEUS-ready input file ( ZEUS.L03), and a series of graphical output files ( CM3DT.L **). 
Of the three, CM3DT.OUT and the graphics files are of interest; ZEUS.L03, although 
critically important in the running of ZEUS, was in an unreadable format to the user. 
The graphic output served as a common sense check of the solution. Using FAST and 
a script written by Edge [Ref. 5], the solution was graphically depicted. As an exam- 
ple, Figure 17 depicts the solution (pressure contours) for the H4 at 0° angle of attack 


at Mach 5. From this graphic output, an inappropriate grid, poorly shaped shock, 





Figure 17. Example of CM3DT Graphic Solution Output 
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or other anomalies such as unexpected pressure changes were easily identifiable. The 
text output file provided many more details and, in essence, walked the user through 
the process. This file was particularly useful for troubleshooting. For example, an 
inconsistency identified in the graphic output could be isolated and examined in the 
test output. 

Finally, ZEUS.LO3 was the critical link between CM3DT and ZEUS. This file 
contained the values of p, p, u, v, and w at each grid point as well as bow shock 
speeds and locations in the z = ZSTART plane. This information is packaged in a 
format that can be read by ZEUS. 
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V. APPLYING THE ZEUS CODE 


A. BACKGROUND 


Having obtained an initial solution plane, the afterbody was analyzed using 
ZEUS. Application of the ZEUS code was greatly facilitated by the use of the ZEUS 
Graphical User Interface (GUI) created by Edge and Clarke [Ref. 5]. The ZEUS GUI 


offered three advantages over normal command line application: 


e easy input file manipulation 


e an ability to observe the computational grid prior to attempting a solution 
run 


e direct links to FAST and pre-written scripts enabling the user to visualize the 
computational grid and the solution 


This section will discuss the application of the ZEUS GUI in obtaining solutions for 
the projectile afterbodies. The methodology, inputs, and output will be presented. 


Except where noted, the discussion in this chapter applies to all cases. 


b. INPUT FILES 
Wardlaw and Davis [Ref. 27] placed all input variables into one file. To 


facilitate the development of the GUI scripts, Edge and Clarke [Ref. 5] split the 


ZEUS input files into six categories: 
e initial conditions 
e integration control 


e zone definitions 


boundary conditions 
e geometry definitions 


e output control 


The specific inputs from each of these categories will be discussed in this section. Àn 


example of these input files can be found in Appendix B. 
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1. Initial Conditions 
The initial condition category describes the freestream conditions used in this 
study. As previously mentioned, ideal gas conditions were applied to this case. The 


inputs used are shown in Table IV. 


Angle of attack | 0-15 
Angle of yaw 


Mach number 
Pressure 
Density 





Table IV. Initial Condition Input Variables 


2. Integration Control 

The integration controls are used in this code to specify the numerical con- 
straints that the code uses. Included in this category are variables describing the 
axial station at which ZEUS ends it calculation in addition to the CFL safety factor 
that determines the marching step size based on a percentage of CFL value. For runs 
where no modifications to the grid or zoning schemes were made, the length of each 
projectile was used for ZETAEND. The step size safety factor value used was the 
default value (9096); the use of this value was based on previous work performed on 
similar objects [Ref. 4, 15]. Decreasing the value of FCFL significantly increased 
the number of steps required to achieve a solution. The maximum number of steps, 
KEND, serves only as a safety in this study. The number of steps required to achieve 
a solution varied based on the complexity of the projectile shape, the angle of attack, 
and the velocity. The H4 at o — 0? and Mach 5 required less than 900 steps. How- 
ever the H4A at a = 15° and Mach 8.2 required over 10,000 steps. In most cases, 
the number of steps required to obtain a flow field solution exceeded maximum value 


allowed in the ZEUS GUI. Continuation of the process beyond the maximum number 
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of steps was performed by restarting ZEUS using the last known solution plane as an 
initial boundary condition, a process identical to that used with the nose solution. 

In Chapter II, application of the limiter, K, in (11.46), was described. The input 
variable for this limiter, XKI, is used. Large values of this variable typically produce 
sharp shocks, but often lead to pressure oscillations. Conversely, small values result 
in less than second order accuracy and induce excessive shock smearing. Wardlaw 
and Priolo [Ref. 27] noted that the default value of 1, produced accurate results and 
that non-default values should be used to “revive” a faulty computation. Finally, 
one input not listed in the Table V, ZSTART, was required by ZEUS; this input, 
like many others, was read by ZEUS from the solution file from CM3DT, ZEUS.L03. 
The value of this variable was specified in CM3DT.INP. The values used are shown 
in Table V. 


ZETAEND | z value for calculation termination | 540 
KEND Maximum steps 


FCFL Step size safety factor 
XKI Limiting constant 
DFAC Turn angle multiplication factor 





Table V. Integration Control Input Variables 


3. Zone Definition Input Variables 

The variables in this input file are one of two files used by ZEUS to construct 
the computational grid. The values chosen for this input file were based on work by 
Priolo and Wardlaw [Ref. 15]. They successfully applied ZEUS to elliptic shaped 
projectiles (2:1, 2.5:1, 3:1 ellipticity ratios) using a single-zoned 72x72 mesh and a 
singled-zoned 36x36 mesh. Use of a denser mesh was chosen to provide greater reso- 
lution particularly on the projectiles with strakes. Because the strakes were modeled 
as part of the body, only a single zone was used. Use of multiple zones in ZEUS are 


typically applied when fins are added to the body. Although use of this feature was 
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originally considered for this problem, the modeling of the projectiles as a body-alone 
object precluded the requirement for multiple zones. À potential use for additional 
zones 1s discussed in the grid sensitivity analysis in Chapter VI. The values and de- 


scriptions of the variables in this file are shown in Table VI. A graphic description of 





Number of zones 1 































NA Number of cells in € direction 12 
MA Number of cells ın n direction [2 
MAZ(IZN) | Number of cells in 7 direction for each zone | 72 


INPANG 





Zone spacing (equidistant=1) 1 


Table VI. Zone Definition Input Variables 


these variables is shown in Figure 18. 
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Figure 18. Zone Definition Variables 
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4. Boundary Conditions Input Variables 

Based on the limitations of ZEUS, default boundary conditions are encoded 
in ZEUS expediting the input setup. Of interest are two input variables, IASYM, 
denoting pitch plane symmetry, and ISHOCK, which defines whether the outer 
boundary is to be shock fitted. Because the projectiles are symmetrical about the 
pitch plane, the pitch plane symmetry boundary condition was used. Based on slender 
body theory, it was assumed that all significant flow information was located between 
the body surface and the bow shock wave. Initial data required to perform the shock 
fitting was obtained from the nose tip solution in ZEUS.L03. 'The information in this 


file overrides an user input values from the GUI. The variables are listed in Table VII. 









IASYM | Pitch plan Symmetry (yes=1) 1 
ISHOCK | Edge 3 to be fitted by calculation of shock (yes=1) | 1 
Table VII. Boundary Condition Input Variables 


5. Geometry Definition Input Variables 

These variables provided the remaining information required to construct the 
computational grid. In addition, information regarding the geometry of the body 
was input through these variables. Of note in the values listed in Table VIII are 
those regarding grid clustering. The grid was clustered about the shoulder of the 
projectile. This was accomplished by doubling the spacing at the outer zone edges (2 
and 4). By doing this, the grid at the shoulder was twice as dense of that at the lee 
and windside meridians. Clustering in the 7 direction was used to most accurately 
predict phenomena that would possibly occur at this location. Priolo and Wardlaw 
[Ref. 15] showed that ZEUS was not able to accurately compute the pressures at the 
shoulder even with a clustered grid. Fortunately, they showed that these inaccuracies 


had little effect on the overall solution. 
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ICORD Type of coordinate system (cylindrical=1) | 1 
IMESHF | Mesh clustering in € direction (yes=1) 
IMESHG | Mesh clustering in 7 direction (yes=1) 
CLUSE4 | Grid spacing at edge 4 boundary 

CLUSE2 | Grid spacing at edge 2 boundary 






















Table VIII. Geometry Definition Input Variables 


6. Output Control Variables 

The majority of the variables in this category were used to format the output 
and are described in [Ref. 5]. However, the reference area variable, AREF, was 
defined in this file. This area was used in the computation of the force and moment 
coefficients. The reference area used in this research was 21,600mm*, the cross- 


sectional area of the tail of the H4. This is the convention used by Edwards and 


Roper [Ref. 10]. 


C. PREPARING THE CM3DT OUTPUT FOR USE IN 
ZEUS 


Having created the required input files, the output from CM3DT had to be 
readied for use with ZEUS. ZEUS.L03, the solution output file from CM3DT, provided 
ZEUS with the flow properties at z = 5mm. However, this solution plane consisted of 
18x24 grid points as defined in CM8DT.INP as NAxMA. The values at these points 
were used to interpolate values for the grid defined by zondim.inp and geo.inp using 
convert.x, part of the ZEUS package and integrated into the GUI. This code took 
the fle RESTART (ZEUS.LO8 was copied to RESTART) and created a new initial 
plane aligned with the mesh characteristics (number of zones and cells, clustering 
,etc.) defined by the user and wrote the output to the file START. Details on the 
interpolation algorithm used in convert.z is available in [Ref. 27]. A simple linear 


interpolation of the previous solution to a new grid is accomplished by re-defining the 
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number of grid points (NA and MA) and establishing new functions describing the 
clustering. The START file is the default file that ZEUS uses when the “Restart” 


option is selected in zntcntrl.inp. 


D. EXAMINING THE COMPUTATIONAL GRID 

The ZEUS GUI offers an extremely powerful tool called Pseudo-ZEUS. This 
program produced the solution grid within minutes (on a Silicon Graphics Octane). 
By producing the solution grid prior to computing a full (and potentially erroneous 
solution), the user was able to make time saving modifications. Since most full so- 
lutions required between 15 to 25 minutes, this program saved significant computer 
time by allowing errors to be corrected prior to the initiation of the full computational 
process. 

Most problems were fixed prior to computation. However, one problem with 
the grid of the models with strakes was not corrected. On one side of each strake, 
orthogonality of the grid is poor. Figure 19 shows the skewness of the grid about this 
strake side. Additionally, only two cells define this side of the strake. Better orthogo- 
nality was achieved by dividing the body into multiple zones. This was performed for 
one case and the results are discussed in the grid sensitivity analysis in Chapter VI. 
An attempt to add more cells resulted in a worse grid structure and an improperly 


defined shock. 


E. ZEUS OUTPUT 

Like CM3DT, ZEUS provided output in visual and numerical form. Typically, 
the visual solution was used for immediate feedback. FAST scripts in the ZEUS GUI 
allowed the user to automatically link into FAST. As an example, Figure 20 shows 
the pressure contours on the H5 at o — 15? at Mach 6. 

The full text output in echoexecut.txt provided the velocity profile and the 


surface pressures. It also showed the ZEUS calculated force and moment coefficients 
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Figure 19. Grid Skewness About the Strake Side 





at each axial station, the last being the total force or moment coefficient. Additionally. 
the center of pressure was provided. The output in echoerecut.trt was the source for 


the analysis in Chapter VI. 


F. A SPECIAL CASE: H4 AND H5 AT o — 157, MACH 7 
AND 8.2 


In each of these four cases. a problem arose in prediction of the shock location 
and therefore the generation of the compuational grid. When the process mentioned 
in the previous sections was applied in these cases, the shock (and grid boundary) was 
severely deformed with an apparent oscillation about the shoulder of the projectiles. 
This problem is shown in Figure 21. 

Investigating the problem revealed that the grid was highly clustered on the 


windward side. Although this grid density caused no problems at lower Mach num- 








Figure 20. ZEUS Graphic Output (H5 at a = 15° at Mach 6) 


bers, it was hypothesized that the combination of high Mach numbers and the density 
of the grid could not be treated accurately by ZEUS’ grid generation algorithm. Al- 
though this claim remains unsubstantiated, the following modifications to the grid 
solved the problem. 

Based on a recommendation by Edge [Ref. 29], a coarser grid (24x24) was 
created and a solution was obtained for the first 20mm of the projectile. At this 
axial station, the grid was returned to its original dimensions and. convert.r was run. 
Pseudo-ZEUS was run for the remaining 520mm showing a properly shaped shock 
and well ordered grid. The result is depicted in Figure 22. 

This problem demonstrates one of the strengths of ZEUS. the ability to stop 
computation before the end of the projectile and redefine the grid (number of cells and 
clustering) and zones. Although this feature was used in this research for a special 


circumstance, its application in the analysis of more complex shapes is obvious. 
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Figure 21. Computational Grid of the H4 at a = 15°, Mach 8.2 (Unmodified) 








Figure 22. Computational Grid of the H4 at a = 15°, Mach 8.2 (Modified) 
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VI. RESULTS 


The goal of this chapter is twofold: 


e to validate the results computed by ZEUS 


e to compare the aerodynamic characteristics of each projectile configuration 


In this chapter, lift, drag, moments, center of pressure, normal force distribution, and 
pressure distribution are analyzed. Particular attention 1s placed on the effects of the 
addition of stabilizing surfaces on these characteristics while generalizing the Mach 
and incidence behavior. À complete listing of the computational results is given in 
Appendix C. 

In this chapter, the non-dimensionalized results are compared. The reference 
area used was 21,600mm?, the planform area, and the reference length used was 
940mm, the projectile length. These are the reference parameters used by Edwards 
and Roper [Ref. 10] and allowed the results from this study to be compared with 


previous results. 


A. VALIDATION 

Validation of the computational results was conducted using experimental and 
computed data published by Edwards and Roper [Ref. 10]. In their work, the H3 
was tested experimentally at Mach 8.2 and over a wide range of angles of attack. 
Additionally, runs were conducted using a parabolized Navier-Stokes (PNS) code, 
CHASM, and Missile DATCOM. For each case, Missile DATCOM was run twice: once 
using second order shock expansion (SOSE) and once using Newtonian corpuscular 
impact theory. The results obtained by ZEUS capture the trends demonstrated by 
the experimental data. For most calculated aerodynamic coefficients, ZEUS tends to 
under-predict their magnitude. However, the results obtained from ZEUS and the 


PNS code generally agree. 


— 
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Figure 23 compares the coefficients of hft. Both CHASM and ZEUS under- 
predict lift whereas Missile DAT COM over-predicts. As angle of attack is increased, 
the difference between the ZEUS and CHASM results becomes more evident (the 
slope of the ZEUS Cy — a curve is smaller than those of other sources). All solvers 
capture the trend depicted by the experimental data with the absolute error for ZEUS 
and CHASM being smaller than those for Missile DATCOM. 
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Figure 23. Comparison of C' versus a Curves 


Figure 24 illustrates that ZEUS significantly under-predicts the drag particu- 
larly at higher angles of attack. This under-prediction can be partially attributed to 
ZEUS' inability to account for boundary layer effects. All solvers under-predict the 
experimental results, however the error incurred by using ZEUS is significantly larger 
than the other solvers. As with lift. the ZEUS’ results are most closely related to the 
CHASM results. Noteworthy is the fact that the absolute error between ZEUS and 


CHASM remains relatively constant throughout the range of incidences. 
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Figure 24. Comparison of Cp versus a Curves 


An analysis of the lift-to-drag ratio curve (Figure 25) shows that ZEUS over- 
predicts the magnitude of the lift-to-drag ratio especially at higher incidences. This 
error, which is partially based on the under-prediction of the drag, remains relatively 
constant throughout the range of incidences. More importantly, it demonstrates 
that the under prediction of drag is more significant than the under prediction of 
lift. Despite its inability to correctly predict the magnitude of the lift-to-drag ratio, 

"n 

Both CHASM and ZEUS are superior at predicting the coefficient of moment 

(Figure 26). Results obtained by ZEUS and CHASM are nearly identical for all 


L 
compared to the other solvers, ZEUS was able to best predict the location of (=) 


incidences. Although both tend to under-predict the moment, the results are in 
better agreement than the Missile DATCOM results which significantly over-predict 


the moment coefficient. 





: L 
Figure 25. Comparison of — versus a Curves 
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В. AERODYNAMIC FORCES AND MOMENTS 


Increasing speed slightly decreases lift for all configurations. Graves [Ref. 
7], obtaining similar results, ascertained that the reduction in lift was caused by 
the diminishing effect of the vortex flow. Increases in incidence result in non-linear 
increases in lift for all projectiles i.e. increasing incidences causes an increase in the 
slope of the Cr, — o curves. This change is the o? variation as noted by Kontis et 
al. [Ref. 9]. Kontis et al. [Ref. 9] also ascertained that the leeside separation of 
the flow seen in schlieren visualization studies caused increased suction on the leeside 
thereby accounting for some of the increased lift. Addition of any of the studied 
stabilizing surfaces results in an increase in lift over the baseline projectile. For all 
Mach numbers and incidences, the flared projectiles, HS and H5A, generated the 


greatest lift. The flared projectiles’ ability to achieve greater lift can be attributed 
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Figure 26. Comparison of Cy versus a Curves 


to a larger planform area. Both Graves and Founier et al. [Ref. 7, 8] noted that 
increased planform area and aspect ratio in elliptical projectiles increased lift. After 
the flared configurations, the H4A had the next largest lift. A study by Kontis et 
al. [Ref. 11] shows that the flow interference caused by the addition of the strakes 
increases lift. This increase in lift can be attributed to an increase in pressure which 
will be discussed in a later section. The difference in lift between the H5 and H5A 
was negligible. A comparison of these configurations shows that the addition of the 
flare appears to dominate. Additionally, because the strakes begin to blend into the 
flare, they create less cross flow interference than that which is found on the H4A. 
Figure 27 graphically depicts these results for a sample case (Mach=8.2). 

As the Mach number increases, there 1s a significant decrease in the drag. In 
supersonic flow, increased speed creates sharper bow shock that are less normal to 


the body surface. As this occurs, the shock weakens and incurs less drag. The effects 
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Figure 27. Cz versus a (Mach 8.2) 


of increasing incidence on drag are similar to those seen on lift. Like the lift results, 
increased incidence causes increased drag. Similarly, Figure 28 shows the effects of 
the addition of control surfaces on drag. A large increase in drag is incurred by the 
addition of the flare. This finding is consistent with the results computed by Guidos 
[Ref. 13] on cone-cylinder-flare configurations. Addition of strakes causes a small 
increase in drag. Larger increases occur at higher incidences based on a greater effect 
on the bow shock [Ref. 11]. 

The relationship between lift and drag is depicted in Figure 29. Increasing 
Mach number results in a significant increase in the lift-to-drag ratio for all configu- 
rations particularly at smaller incidences. Generally, (L/D)mar occurred between 5° 
and 10° incidence. (L/D)mar was achieved at lower incidences as the Mach number 
was increased. At non-zero angles of attack, the projectiles with flared tails have 


significantly lower lift-to-drag ratio indicating that the increase in drag incurred by 
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Figure 28. Cp versus o (Mach 8.2) 


adding the flare is much greater than the added lift obtained. The addition of strakes 
on the H4 causes little change in the lift-to-drag ratio. This is inconsistent with the 
findings of Kontis et al. [Ref. 9] who found that the addition of strakes resulted in a 
17% increase in the lift-to-drag ratio. Again, this can be attributed to ZEUS’ under 
prediction of drag. 

Generally, the pitching moment decreases with increasing Mach number and 
increases with increased incidence for all configurations. Figure 30 shows the effects 
of the stabilizing surfaces. The addition of any studied surface increases the moment. 
The flare significantly increases the moment, but adding strakes to the projectile 
with a flare tail has little impact on the moment. The addition of the strakes alone 
does increase the moment compared to the body alone configuration. However, this 
increase is not nearly as dramatic as that found by adding a flare. As expected, the 


configurations demonstrating the highest stability incur the greatest drag. 
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Figure 29. = versus o (Mach 8.2) 


C. PRESSURE DISTRIBUTION AND CENTER OF PRES- 
SURE LOCATION 


Surface pressures were calculated for each grid point. The following section 
analyzes the surface pressures on the windside (¢ = 0°) and leeside (¢ = 180°) 
meridians. 

For all configurations, increasing angle of attack causes a decrease in the leeside 
pressure and an increase in the windside pressure. Increasing Mach number has a 
similar effect. Figures 31, 32, and 33 compare the wind and leeside pressures of the 
H5, H4A, and H5A, respectively, to those of the H4. These comparisons are over the 
range of incidences at Mach 8.2. 

At all studied angles of attack, there is a significant increase in both the 
windside and leeside pressures when the flare is added. The pressures acting on the 


H4 and H5 are identical up to the vicinity of the body-flare junction. The largest 
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Figure 30. Cy versus o (Mach 8.2) 


absolute increases in pressure occur at higher angles of attack. However, in relative 
terms, the pressure increases at lower angles of attack are larger. It is important 
to note that Hughson and Chapman [Ref. 12] showed that the magnitude of the 
pressures just before the flare-body junction are significantly different depending on 
the nose shape; blunt nosed configurations caused pressures at this location to be 
dramatically larger. 

Application of strakes caused an increase in windside pressures and a decrease 
in the leeside pressures. These pressure changes are moderate compared to the rel- 
atively large increases seen in the flare application. Additionally, they appear to be 
less dramatic than the pressure differences seen in the application of strakes to a right 
elliptic cone Kontis et al. [Ref. 9]. As expected, increasing angle of attack results in 
pressure increases and decreases on the wind and leesides respectively. 


Use of both strakes and a flare result in a combination of the previously men- 
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Figure 31. Leeside and Windside Pressure Distributions of the H4 and H5 at Mach 
D 

tioned effects. Again a dramatic pressure increase is seen at the body-flare junction. 
Increasing angle of attack had the same effect as in the other configurations. Figure 
34 compares the pressures on the H5 and H5A to emphasize the effects of the strakes 
on the flare. Like the effect on the aerodynamic forces and moments, the flare is 
dominant in the effect of surface pressure. The application of strakes has little effect 
on the pressure distribution particularly on the windside. The center of pressure was 
measured as a percentage of the total projectile length. Table EX shows the effects 
of the control surfaces on the center of pressure. The addition of the flare had the 


largest impact on the center of pressure resulting in a nearly 5% change. The addition 
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Figure 32. Leeside and Windside Pressure Distributions of the H4 and H4A at Mach 
8.2 

of strakes had little effect on the center of pressure. These results do not agree with 
those of Kontis et al. [Ref. 9] who found that the addition of strakes caused a 4% 


change in the shift (to the rear) in the center of pressure with the addition of strakes. 


D. NORMAL FORCE DISTRIBUTION 


The following section analyzes the effect of each projectile geometry on the 
normal force distribution. 


Generally. increasing Mach number results in a decrease in the normal force. 
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Figure 33. Leeside and Windside Pressure Distributions of the H4 and H5A at Mach 
8.2 


Center of Pressure (% of Projectile Length) 
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Table IX. Centers of Pressure (Mach 8.2) 
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Figure 34. Leeside and Windside Pressure Distributions of the H5 and H5A at Mach 
5.2 

This decrease in normal force is more dramatic at higher angles of attack. Addition- 
ally, higher angles of attack increase the normal force. However, neither speed nor 
incidence change the normal force distribution trend. 

The projectile geometry did have an effect on the normal force distribution. 
Figure 35 depicts the normal force distribution on each of the projectiles at Mach 8.2 
and @ = 5°. The general shapes of the force distributions are nearly identical for all 
configuarations for the first 430mm, the location of the body-flare junction. As ex- 
pected, the straked configurations generated greater normal force, yet the distribution 


trends were similar to those without strakes. 
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Figure 35. Normal Force Distributions at Mach 8.2, a = 5° 


At the body-flare junction, a dramatic increase in normal force acting on 
the projectiles with flares occurs. Figure 36 shows the normal force distribution for 
z = 430mm to z = 540mm. As expected, as the flares taper into the tail of the H5A, 
the difference in the magnitude of the normal forces become smaller. Also noteworthy 
is the fact that the normal forces acting on the H4A are larger than those acting on 


the H5 until z 2 470mm, over a third of the length of the flare. 


E. PITCHING MOMENT DISTRIBUTION 


The pitching moment distribution was analyzed for each projectile configura- 
tion. Figure 37 compares the pitching moment distribution for each configuration 
at Mach 8.2 and a = 5°. Like the normal forces, the pitching moments are nearly 
identical for all configurations from the nose to the vicinity of z = 250mm. At this 
axial station, the effects of the strakes become apparent causing a moderate increase 
in the moment. Moving aft. the difference between the H4A’s and the H4's moment 
increases demonstrating the influence of increasing strake height on the moment. The 


moment distribution of the H4 and H5 are identical until the vicinitv of the flare-body 
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Figure 36. Normal Force Distributions on the Tail Section at Mach 8.2, a = 5° 
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Figure 37. Pitching Moment Distributions at Mach 8.2, a = 5° 
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junction. Comparing the H4A and H5A shows similar results. At the flare-body junc- 
tion, the moment of the flare-tailed projectiles increases significantly. The H5A has a 
greater moment than the H5 until roughly z = 510mm when the effects of the strakes 
become negligible. Similar to the results of the normal force distribution analysis, 
the pitching moments acting on the H4A are larger than those acting on the H5 until 


= — 4 (mm! 


F. GRID SENSITIVITY ANALYSIS 

As mentioned in Chapter V, problems arose in the grid generation for the 
configurations with strakes; the grid around the strakes was skewed. Based on the 
preliminary nature of this study, this grid was used for all runs. However, an improved 
grid was constructed and tested on one case. This section will describe the results of 
this analysis. 

The construction of the improved grid was more detailed than the single-zoned, 
clustered mesh used throughout this research. The improved grid was comprised of 
five zones: one covering each strake and three covering the body surface. The ZEUS 
grid generation algorithm was unable to properly capture the bow shock shape so the 
mesh was made uniform and more coarse (35x35). 

Figure 38 shows the improved grid around the strake on the H4A. A test run 
was conducted on the H4A at Mach 8.2, a = 5°. The results are shown in Table X. 
The analysis shows that the improved grid had almost no impact on the calculation 
of the axial force and a small effect on the normal force. The improved grid resulted 
in a roughly 2% increase in the normal force. Of the coefficients calculated by ZEUS, 
the one most affected by the improved grid was the pitching moment coefficient. 
Use of the new mesh increased the pitching moment coefficient over 2%. The center 
of pressure shifted slightly toward the base with the improved grid. Overall, the 
improved grid had a small impact on the results. 


A second grid sensitivity analysis was conducted to attempt to account for the 
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F igure 38. Improved Grid for the H4A 


under prediction of drag. Because the mesh is “stretched” on the leeside at higher 
angles of attack, the mesh was clustered normally. The analysis was conducted on the 
H4 at all studied angles of attack resulting in negligible differences with the previously 


computed results. 


63721-03 | 1.8578Е-02 | 5.5942Е--00 
Improved | 3.3298E-02 | 2.6288E-03 | 1.9095E-02 | 5.9665E-+-00 
TASK -0.32% ZA 1.22% 


Original | 3.2571Е-02 | 2. 
2 





Error 
Table X. Grid Sensitivity Results 
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VII. CONCLUSIONS 


Four projectile configurations with elliptical cross-sections have been analyzed 
using a spatial marching Euler solver. The results were compared to available exper- 


imental and computed data. The evaluation of the results shows: 


e ZEUS is able to reasonably predict trends in aerodynamic characteristics for 
these projectiles. Neglected boundary layer effects and difficulty constructing 
the computational grid about the strakes cause error often leading to poor 
agreement with experimental results. ZEUS is capable of predicting lift and 
moment to a level of accuracy usually found in the more computationally 
intensive PNS codes. 


e Generally, the effects caused by changes in Mach number and moderate changes 
in angle of attack are similar for all studied configurations. 


ə The application of strakes, a flare, or both cause an increase in lift. The 
current flare design induces greater lift than the strakes. 


e Addition of a flare significantly increases drag. Although the strakes also 
generate additional drag, the amount of drag created by the flare is much 
greater. 


e For all configurations, the maximum lift-to-drag ratio occurs at approximately 
5°. The drag penalty incurred by increasing lift is largest for the projectiles 
with flares. The difference in lift-to-drag ratio for designs with and without 
strakes is very small. 


e The configurations with the flare had the highest pitching moment coefh- 
cient. Adding strakes increases the coefficient of pitching moment regardless 
of whether a flare is present. 


e Addition of strakes causes increased pressure on the windside and a decrease in 
pressure on the leeside. Application of a tail flare results in a dramatic increase 
in pressure on both the wind and lee sides from the flare-body junction to the 
base. 


e Normal force distribution trends are identical for all projectiles until the body- 
flare junction. The flare causes a large increase in the force distribution. This 
Increase is more dramatic for the flared projectile without strakes. 
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e The pitching moment for the flared configurations increases significantly start- 
ing at the body-flare junction. As the strakes become more prominent on the 
projectiles (toward the tail of the H4A or the body-flare junction of the H5A), 
the pitching moment increases. 
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APPENDIX A. CODES USED FOR THE 
CONSTRUCTION OF THE PROJECTILE 
GEOMETRY FILES 


The following are MATLAB codes designed to create files describing the geom- 
etry of each projectile. The output of each code is a file in PLOT3D format that can 
be used in CM3DT and ZEUS. The codes shown create half of the projectile (to be 
used for incidence computations only). Additionally, the nose generation codes have 
been omitted. But for several modifications the nose generation codes are nearly 
identical to the projectile generation codes. As noted in Chapter III, the following 


modifications were made to full projectile codes: 


e The projectile length is reduced to 10mm 


e The number of points defining to nose are reduced in the circumferential di- 
rection and increased in the axial direction 


e The results are output to coincide with the CM3DT coordinate system 


1. THE H4 GENERATOR-Z4GENERATOR.M 
/h4generator.m 

AThis code creates a .P3D file for an elliptical projectile. 
A3MAY2000 


даата-П;Х-П;Ү-П;2-П; 


z-linspace(0,540,109); 
theta=linspace (pi/2,-pi/2,50); 
for 3=1:50 
for i=1:109; 
a-sqrt(15.*z(i))/3; 
besgrt(159*z(0320)75. 
y-a*cos (theta(j)) ; 
x-b*sin(theta(j)); 
X=[X x];Y=[Y y]:Z2=[Z 2G:)] 
end 


о 


end 


one-1; 
jmax=j; 
imax-109; 
data=[X ,Y,Z]; 


fid = fopen(’h4.grd’,’w’); 

fprintt fid, fd ове) 

fprintiltid, 4d жача хп лах ушак опе): 
fprintf(fid, 412.8e /1295e 2 12 s o и атау 
fclose(fid); 


2. THE H4A GENERATOR-HZAGENERATOR.M 


),h4agenerator.m 
“This code creates a .P3D file for an elliptical projectile. 
%3MAY2000 


data-[] ;X12[] ; Y12 [] ;Z12 [] ; X2 [1]; Y2 []; Z2 [] ; imax20; |шах-0; 
z-linspace(0,540,109); 


hI. Creates body shell from bottom(pi/2) to the bottom, right strake 
theta=linspace(pi/2,.9614,10); 
for 1=1:10 
a-sqrt(15.*z)/3; 
b=sqrt (15.*z)/5; 
y=axcos(theta(i)); 
x=b*sin(theta(i)); 
X=(X,xJ ;Y=(Y,y];Z=(Z,z]; 
end 
imax-imax*i; 


AII. Creates bottom,right strake 

АА. Creates side from surface to top of strake (straight line). 
%This is accomplished by determining the equation of the line at each 
Zvalue of z. This line is broken into 4 intervals (size h) vane 

„еһе point for (y=y1+ih,x=my+b) is determined for all z. 

phi=.8844; % This the angle to the strake top point. 

al=sqrt (2.268. *z) ; 

bi-sqrt(.9833.*z); 

yi=ai.*cos(phi) ; 


x1=b1.*sin(phi); 


forsr-1:3 
for j=1:109 
if j==1%this conditional statement eliminates Singularity at 
КОТ ырп 
х=0;у=0; 
x= Dax] ;Y=[Y,y] 22200]: 
else 
m=(x1(j)-X(j+981))/Gy10 ESE 
bint=x1(j)-m.*y1(j); 
h=abs (y1(j)-Y(j+981))/4; 
y=Y(j+981)+(h*i); 
x=m*y+bint ; 
X=[X,x] ;Y=[Y,y] ;Z=[Z,z(j)]; 
end 
end 
end 
imax-imax-*1i; 


АВ. This loop will construct ohe top of the strake. 
phi=linspace(phi, .7158,5); 
for 1=1:5 
y=al.*cos(phi(i)); 
x=bi*sin(phi(i)); 
X=[X,x] ;Y=[Y,y] ;Z=[Z,z] ; 
end 
lmax=imaxti; 


АС. This loop will construct the other side of the strake. This is 
Zan modeling approximation of the actual shape (straight). The shape 
Ais modeled with a circular arc varying arc length. The circle from 
Awhich this is derived has a varying center and radius that are 
4calculated in the code. 

eta=.5203;%angle to dummy (modeled) point on body surface 
zeta-.6578;Áangle to actual point on body surface 

y2=a.*cos(eta); 

x2=b.*sin(eta); 

y3za.*cos(zeta); 

x3=b.*sin(zeta); 


for 1212.10 
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forn) —1 2109 
if j==1/this conditional statement eliminates singularity at 
horigin 
X7 [X,0] ; Y [Y,0] ; Z2 [2,0] ; 
else 
AApproximate three points of circle 
ymidbad=(y2(j)-Y(j+1853))/2+Y(j+1853) ; 
xmidbad=(x2(j)-X(j+1853))/2+X(j+1853); 
ymidba=(y3(j) -Y(j+1853) ) /2+Y(j +1853) ; 
xmidba= (x3 (j) -X(j+1853) )/2+X (5 +1853) ; 
ythird=(ymidbad-ymidba) /2+ymidba; 
xthird=(xmidbad-xmidba) /2+xmidba; 
“Calculate center and radius of circle 
x4=X(j+1853); 
у4=Ү(ј+1853); 
h2(C(xthird^2-x2(j)^2)*(ythird^2-y2(j)^20)/... 
а )) (G20( о 
(y2()2) 2-y2"2))/ 2*(X2CGJ)=x4)))/ а уа ee 
CASAL AA 
КЕСИ 2 AD A Y ee 
(2*(x2(j)-x4))+((y4-y2(j))/(x2(j)-x4))*h; 
pa ERAS AI SA h)°2)5 
deltax=abs (x2(j)-x4)/10; 
x=x4+deltax*1; 
Vo aoqre (lr 20x kl 2) eh: 
X-[X,x] ; Y2[Y, y1;Z* [Z,z C2]; 
end 
end 
end 
lmax-imax-i; 


AIII. Creates projectile shell between strakes 
theta-linspace(eta,-eta,20); 
for 1=1:20 
yza*cos(theta(i)); 
x-b*sin(theta(i)); 
Xo [Xj xl ҚА УА ҮЛҮЛЕ 
end 
lmax-imax-*i; 


AIV. Creates top strake. Since this strake is a reflection of 
Athe previously modeled strake reflected about the y-axis. 
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7921-7 28 

j=3052- (i*109) +1; 

XML 216108)] ;Y=[Y, 76 Oey Z= (2207-74108) 9; 
end 
imax=1imax+1; 


one=1; 

data=[X,Y,Z]; 

jmax=109; 

fid = fopen(’h4al.grd’,’w’); 

Bprintf(fid, Ad An? one); 

Dprintf(fid, Ad Xd Ad Nn) jmax im as one), 
torintf(fid ^ 419-8e 419 98e MT? 8e Nn? data), 
fRellose (tid) 


3. | THE H5 GENERATOR-Z5GENERATOR.M 
Ah5generator.m 

AThis code creates the H5 elliptical projectile with flared tail. 
403MAY 2000 


data-[] ;X2[] ; Y2[] ;Z- [.] ; 
theta-linspace(pi/2,-pi/2,50); 
bigx-35*cos(theta); 
bigy=22.85*sin(theta); 


tor y=1:50 

for z=0:5:430 
a-sqrt(15.*z)/3; 
b=sqrt(15.*z)/5; 
x-a*cos(theta(j)); 
y=b*sin(theta(j)); 
Xs [X,x]; Y2[Y, y1;Z- [Z,z] ; 

end 


dxdz-(bigx(j)-x)/(540-430); 
dydz=(bigy(j)-y)/(540-430); 
bxz=x-dxdzx*430; 
byz-y-dydz*430; 


for 1=435:5:540 
x=dxdz*i+tbxz; 
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y*dydz*itbyz; 
Z=1; 
X=[X,x];Y=[Y,y];Z=[Z,z); 
end 
end 


one=1; 
jmax=j ; 
іпах-109; 
data=[Y IX 7]; 


fid = Topen( hS gerd uF 

fprintí fid, 7a none). 

fprintrierid.z2. rede dsa nn ymax, max one); 

Tb se 12 88 Ута ве data 
fclose rior 


4. THE H5A GENERATOR-H5AGENERATOR.M 
hh5agenerator .nm 

‘This code creates the H5 elliptical projectile with flared tail. 
LO3MAY 2000 


data=[1;X=[] ;Y=[] ;Z=[] ; imax=0; jmax=0 ;x2=[] ;y2=[] ;x3=[] ;y3=[]; 


theta-linspace(pi/2,.9614,10); 
bigy=35*cos(theta); 
bigx=22.85*sin(theta); 


for 0 

for z=0:5:430 
a=sqrt(15.*z)/3; 
b=sqrt(15.*z)/5; 
y=a*cos (theta(i)); 
x=b*sin(theta(i)); 
X= (Xe x Zz=lz 215 

end 


dxdz=(bigx(i)-x)/ (540-430) ; 
dydz=(bigy(i)-y)/ (540-430) ; 
bxzZz=x-dxdz*430 ; 
byz-y-dydz*430; 


for j2435:5:540 
x=dxdz*j+bxz ; 
y=dydz*j+byz; 
251; 
Хә к x]; YS [Yo yd ZS EZ. 2]; 
end 
end 
1max=1maxti; 


z-linspace(0,540,109); 

phi=.8844; ^ This the angle to the strake top point. 
ai-sqrt(2.268.*z); 

bi-sqrt(.9833.*z); 

yl=al.*cos(phi); 

x1=b1.*sin(phi); 


for 1=1:3 
for j=1:109 
if j==1%this conditional statement eliminates singularity at 
horigin 
х=0;у=0; 
X=[X,x];Y=[Y,yJ];Z=[Z,z(1)]; 
else 
m=(x1(j)-X(j+981))/(y1(j)-Y(j+981)); 
bint=x1(j)-m.*y1(j); 
h=abs (y1(j)-Y(j+981))/4; 
y=Y(j+981)+(h*i); 
x=m*y+bint; 
X=[X,x] Y=[Y yJ]:Z=[2Z 2(j)); 
end 
end 
end 
imax-imax-*ti; 


AB. This loop will construct the top of the strake. 
phi=linspace(phi, .7158,5); 
for i=1:5 
y=al.*cos(phi(i)); 
x-bi*sin(phi(i)); 
X-[X, x]; Y2 [Y, y]1;Z- [Z, z] ; 
end 
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imax-imax-*i; 


AC. This loop will construct the other side of the strake. This 
Ais an modeling approximation of the actual shape (straight). The 
ÁAshape is modeled with a circular arc varying arc length. The 
heircle from which this is derived has a varying center and radius 
Athat are calculated in the code. 

eta-.5203;Áangle to dummy (modeled) point on body surface 
zeta=.6578;/angle to actual point on body surface 


for z=0:5:430 

a=sqrt(15*z)/3; 

b=sqrt(15*z)/5; 

y=a*cos(eta); 

x=b*sin(eta); 

yy=a*cos (zeta) ; 

xx=b*sin(zeta) ; 

x2=[x2,x] ;y2=[y2,y] ;x3=[x3,xx] ;y3=[y3,yy]; 
end 


bigy-35*cos(eta); 
bigx-22.85*sin(eta); 
bigyy=35*c0s (zeta); 
bigxx=22.85*sin(zeta); 
dxdz=(bigx-x)/(540-430); 
dydz=(bigy-y)/(540-430) ; 
bxz=x-dxdz*430; 
byz=y-dydz*430; 
dxxdz=(bigxx-xx) / (540-430); 
dyydz=(bigyy-yy) / (540-430) ; 
bxxz=xx-dxxdz*430; 
byyz-yy-dyydz*430; 


for j=435:5:540 

x=dxdz*j+bxz; 

у=ауа2 92 

xx=dxxdz*j+bxxz; 

yy-dyydz*)*byyz; 

x22[x2,x] ; y22 [y2, y]; x32 [x3,xx] ; y32 [y3, yy] ; 
end 


z-linspace(0,540,109); 
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for i=1:10 
for j=1:109 
if j==1%this conditional statement eliminates singularity 
hat origin 
X [X .0];Y=[Y 0] 2=|[2 0]. 
else 
/Approximate three points of circle 
ymidbad=(y2(3)-Y(5+1853))/2+Y(3+1853); 
xmidbad=(x2(j)-X(j #1653) )724k() +1853) ; 
ymidba=(y3(j)-Y(j+1853) ) /2+Y (j +1853) ; 
xmidba=(x3(j)-X(j+1853) ) /2+X (j +1853) ; 
ythird=(ymidbad-ymidba)/2+ymidba; 
xthird=(xmidbad-xmidba)/2+xmidba; 
ACalculate center and radius of circle 
х4-Х(1%1853); 
у4-Ү(1%1853); 
h-(((xthird^2-x2(j)^2)*(ythird^2-y2(j)^2))/... 
(2*(xthird-x2(j)))-((x2(j) "2-x4°2)+(y2(j)"2-y4°2))/... 
(2*(x2.Gj )=x4) ))/ CG 4-4 2G) 7 (x20) exa» thurds e 
ОЙО Еа БР E 
k=((x2(j) "2-x4°2)+(y2(j)"2-y4°2))/(2*(x2(j)-x4))*... 
((y4-y2 (322 / (x2 (j) -x4)) *h; 
r=sqrt((x2(J))-kK) 2202 0 ЕЕ 
deltax=(x2(j)-x4)/10; 
x=x4+deltax*1; 
Еа Ек CE 
if j<109 
Xx=[X,x];Y=[Y,y] ;Zz=[Z,z(j)]; 
else 
theta=asin(x/22.85); 
X=[X,x];Y=[Y,35*cos(theta)];Z=[Z,z(3)]; 
end 
end 
end 
end 
imax-imax-*i; 


AIII. Creates projectile shell between strakes 
theta-linspace(eta,0,10); 

bigy-35*cos(theta); 

bigx-22.85*sin(theta); 


for 1=1:10 

for z=0:5:430 
a=sqrt(15.*z)/3; 
b=sqrt(15.*z)/5; 
y=a*cos(theta(i)); 
x-b*sin(theta(i)); 
у= Е y] z- [2 25 

end 


dxdz=(bigx(i)-x)/(540-430); 
dydz=(bigy(i)-y)/(540-430); 
bxz=x-dxdz*430; 
byz=y-dydz*430; 


for j=435:5:540 
x=dxdz*3+bxz ; 
y=dydz*j+byz; 
2=); 
X=[X,x]; Y=[Y,y];2=[Z,z]; 
end 
end 
imax-imaxt*i; 


AIV. Creates top strake. Since this strake is a reflection 
hof the previously modeled strake reflected about the y-axis. 
for i=l. 3/ 
k=4034-(i*109); : 
X=[X,-X(k:k+108)];Y=[Y,Y(k:k+108)];Z=[Z,Z(k:k+108)]; 
end 
for p=1:109 
Y(p)=0; 
Y(p+8064)=0; 
end 
imaxzimax-*i; 


one=1; 
jmax=109; 
datasz[X,Y,Z]; 


Bud - ғореп /һоа рта” Тұла 


fprintf(fid,*fd Ха” ,опе); 
вре Сеа Ла ха Ха” max max one); 
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fprinti fid 412.8e /12.8e minio ni, data); 
fc ес:а): 
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APPENDIX B. SAMPLE CM3DT AND ZEUS 
INPUT FILES (H4 o ^ 0" MACH 5) 


1. CM3DT INPUT FILES 


a. Geometry Pre-Processor Input File-preproc.inp 
$input 
isa=1,lmax=9,fname=’h4anose.inp’, 
$end 


b. CM3DT Solver Input File~-CM3DT INP 
CM3DT 
$INPUT 
AMINF=5, 
ILAM=1, 
NMAX=11,MMAX=18,LMAX=9, 
KMAX=1000 ,KOUT=1000, 
М5ЕС-0,ХВ0-0.0, 
kense(4)=1, 
rn=.5,xbb=9*8 ,alpha=0.0, 
ktape=1000,kense(5)=1 
$END 
$ZEUS 
na=18 ,ma=24,zstart=5.00,z0=0.0,ibl=1, 
x0=0.0 
$END 


2. ZEUS INPUT FILES 


a. Initial Conditions Input File—initc.inp 
.0 ALPHA 
.0 BETA 
.0 XMINF 
.0 PINF 
.00001 DINF 


© к л Q' € 


b. Integration Control Input File-intentrl.inp 
1 IMOD 
5.00 ZSTART 
540.00 ZETAEND 
1500 KEND 
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SQOMISCRI. 
.00 XKI 
IAPR 
.90 DFAC 
IVIS 


O Orere O 


c. Zone Dimensions Input File—zondim.inp 
1 IZN 
72 NA 
72 MA 
72 MAZ(IZN) 
O INPANG 
ANGMAZ (IZN) 


d. Boundary Conditions Input File—bdc.inp 
1 IASYM 
O NXKE 
O NSUR 
1 ISHOCK 
50.00 1.25 540.75 RSHOCK PHI1SH PHI2SH 


e. Output Control Input File—outentrl.inp 
1 IEPORCE 
1696.46 AEREF 
ЕЕЕ 
1 IPIETSD 
20 INCP3D 
1 “ІШЖЕУӘ 
500 TERINT 
Jo NSRIP 
Tos ES 
SU RE 
加 ER 
50 KSKIP 
100 NPRT 
100 DELZA 
1000 IPLOT 
200 DELZC 
0 ІРІОТМ 
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9.60171E-06 
3.14319 E-02 
8.97250E-02 
1.57475E-01 


4.03579E-05 
3.23951E-02 
8.48458E-02 
1.49157 E-01 


6.02926E-05 
3.18681E-02 
8.04736 E-02 
1.42990E-01 


6.91120E-05 
3.05624 E-02 
7.63588E-02 
1.37767E-01 


9.60171E-06 
3.08883E-02 
8.74438E-02 
1.50349E-01 


4.03579E-05 
3.19527E-02 
8.27962E-02 
1.42506E-01 


6.02926E-05 
3.14881E-02 
7.85803E-02 
1.36660E-01 


6.91120E-05 
3.02309E-02 
7.45924E-02 
1.31696E-01 


ZEUS 


5.31056E-03 
4.86514E-03 
5.28713E-03 
6.80190E-03 


3.91452E-03 
3.66156E-03 
4.38042E-03 
6.05958E-03 


3.07053E-03 
2.96838E-03 
3.86242E-03 
5.63259E-03 


5.26036E +01 
5.33295Е--01 
7.54232E+01 
1.14896E+02 


lil 


5.31056E-03 
7.58610E-03 
2.07874E-02 
4.732717 E-02 


3.91452E-03 
6.47105 E-03 
1.90472E-02 
4.44578E-02 


3.07053E-03 
5.73457 E-03 
1.77778E-02 
4.24492E-02 


2.43535E-03 
5.12325E-03 
1.66983E-02 
4.07946E-02 


2.43380E-05 
1.70176E-02 
5.11541E-02 
9.11934E-02 


4.59032E-05 
0.018046296 
0.048651492 
0.086538066 


7.07231E-05 
0.018026492 
0.046272805 
0.083029835 


9.6953E-05 
0.017462706 
0.043980967 
0.080021005 


APPENDIX C. COMPUTED RESULTS FROM 


1.36877E+03 
2.92363E+02 
3.07865E+02 
3.12715E+02 


6.14199E+02 
3.00817E+02 
3.09642E+02 
3.13297E+02 


6.33418E+02 
3.05457E+02 


3.10502E+02 
3.13562E+02 


7.57533E+02 
3.08545E+02 
3.11028E+02 
3.13655E+02 





7.80023E-05 
3.53338E-02 
9.99963E-02 
1.74542E-01 


1.39519E-04 
3.68539E-02 
9.58653E-02 
1.67111E-01 


1.90365 E-04 
3.68175E-02 
9.21259E-02 
1.61716E-01 


7.80023E-05 
3.45212E-02 
9.70130E-02 
1.65774E-01 


1.39519E-04 
3.62023E-02 
9.31848E-02 
1.58873E-01 


1.90365E-04 
3.62621E-02 
8.96366E-02 
1.53814E-01 


8.45546E-03 
7.78014E-03 
8.43176E-03 
1.08967E-02 


6.28116E-03 
5.86782E-03 
7.04940E-03 
9.82764E-03 


4.95074E-03 
4.76458E-03 
6.27528E-03 
9.23907E-03 


8.45546E-03 
1.08301E-02 
2.56678E-02 
5.57002E-02 


6.28116E-03 
9.05753E-03 
2.35891E-02 
5.27443E-02 


4.95074E-03 
7.95531E-03 
2.21774E-02 
5.07794E-02 


8.25882E-05 
2.05997E-02 
6.05754E-02 
1.06854E-01 


1.31353E-04 
2.21267 E-02 
5.87419E-02 
1.02990E-01 


1.84084E-04 
2.25495E-02 
5.69305E-02 
1.00173E-01 


5.71749E+02 
3.14823E+02 
3.27119E+02 
3.30584E+02 


5.08393E+02 
3.24210E+02 
3.30888E+02 
3.32799E+02 


5.22185E+02 
3.30731E+02 
3.33700E+02 


3.34499E+02 


0 
5 
10 
15 


2.35248Е-04 
3.60126E-02 
8.86731E-02 
1.57282E-01 


2.35248E-04 
3.55291E-02 
8.63288E-02 
1.49637 E-01 


3.9359 7E-03 
3.9754 7E-03 
5.74273E-03 
8.83181 E-03 


3.93597E-03 
7.09905E-03 
2.10534E-02 
4.92384E-02 


2.43470E-04 
2.24378E-02 
5.52365E-02 
9.78781E-02 


5.58872E+02 
3.36448E+02 
3.36377E+02 
3.36046E+02 





1.33857 E-05 
3.61795 E-02 
9.45856E-02 
1.65864E-01 


6.43222E-06 
3.55308E-02 
8.97375 E-02 
1.57426E-01 


3.10100E-06 
3.43067 E-02 
8.53019E-02 
1.51519E-01 


1.42867 E-06 
3.25705 E-02 
8.12259E-02 
1.46669E-01 


1.33857 E-05 
3.55745 E-02 
9.21578E-02 
1.58380E-01 


6.43222E-06 
3.50472E-02 
8.75606E-02 
1.50452E-01 


5.63023E-03 
5.36176E-03 
5.70616E-03 
7.08005E-03 


4.11024E-03 
3.99755E-03 
4.68514E-03 
6.22171E-03 


5.63023E-03 
8.49461E-03 
2.20441 E-02 
4.97677E-02 


4.11024E-03 
7.07906E-03 
2.01967 E-02 
4.67547E-02 


3.10100E-06 | 3.19044E-03 | 3.19044E-03 


3.38963E-02 
8.32954E-02 
1.44880E-01 


1.42867E-06 
3.22167E-02 
7.93569E-02 
1.40293E-01 


3.21061 E-03 
4.09174E-03 
5.70671E-03 


6.18841 E-03 
1.88421E-02 
4.47284E-02 


Mach 8.2 


2.49925E-03 
2.63716E-03 
3.65685E-03 
5.32412E-03 
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2.49925E-03 
5.46583E-03 
1.77060E-02 
4.31034E-02 


1.63766E-07 
1.97623E-02 
5.38735 E-02 
9.57605 E-02 


3.31669E-07 
1.98006E-02 
5.14367 E-02 
9.10322E-02 


9.52992E-08 
1.93819E-02 
4.89908E-02 
8.77040E-02 


1.54847E-07 
1.85784E-02 
4.67130E-02 
8.49345E-02 


6.60657E+00 
2.94963E+02 
3.07570E+02 
3.11766E+02 


2.78444E+01 
3.00931E+02 
3.09523E+02 
3.12257E+02 


1.65951E+01 
3.05078E+02 
3.10135E+02 
3.12568E+02 


5.85283E+01 
3.08019E+02 


3.10554E+02 
3.12708E+02 





3.93450E-02 
1.01158E-01 
1.75503E-01 


7.20787E-06 
3.89064E-02 
9.68843E-02 
1.67941E-01 


4.64259E-06 
3.79903E-02 
9.29519E-02 
1.62732E-01 


2.76933E-06 
3.65974E-02 
8.93852E-02 
1.58575E-01 


3.85136E-02 
9.81789E-02 
1.66824E-01 


7.20787E-06 
3.82469E-02 
9.42146E-02 
1.59819E-01 


4.64259E-06 
3.74329 E-02 
9.04836E-02 
1.54963E-01 


2.76933E-06 
3.61170E-02 
8.70725E-02 
1.51073E-01 


8.19352E-03 
7.82106E-03 
8.30778E-03 
1.04264E-02 


6.00361E-03 
5.86806E-03 
6.89759E-03 
9.26991E-03 


4.6897 2E-03 
4.73671E-03 
6.08162E-03 
8.59352E-03 


3.69011E-03 
3.91411E-03 
5.49806E-03 
8.10824E-03 
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8.19352E-03 
1.12204E-02 
1.81834E-02 
5.54946E-02 


6.00361E-03 
9.23664E-03 
1.63611E-02 
5.24203E-02 


4.68972E-03 
8.02976E-03 
1.51636E-02 
5.04188E-02 


3.69011E-03 
7.08889E-03 
1.42332E-02 
4.88742E-02 


1.74144E-06 
2.24772E-02 
5.97174E-02 
1.04414 E-01 


1.14270E-06 
2.27715E-02 
5.78242E-02 
1.00491E-01 


1.48262E-06 
2.26586E-02 
5.58558E-02 
9.78095 E-02 


1.59348E-06 
2.21796E-02 
5.40592E-02 
9.56867E-02 


6.13806E+01 
3.08493E+02 
3.18782E+02 
3.21271E+02 


8.56085E+01 
3.16056E+02 
3.22292E+02 
3.23121E+02 


1.72450E+02 
3.22073E+02 
3.24492E+02 
3.24564E+02 


3.10717E+02 
3.27263E+02 
3.26587E+02 
3.25844E+02 
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